File indexing completed on 2024-04-06 12:09:06
0001
0002
0003
0004
0005
0006 #include "CommonTools/TriggerUtils/interface/GenericTriggerEventFlag.h"
0007 #include "DQM/TrackingMonitor/interface/dEdxAnalyzer.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "MagneticField/Engine/interface/MagneticField.h"
0012 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0013 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0014 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0015 #include <string>
0016 #include "TMath.h"
0017
0018 dEdxAnalyzer::dEdxAnalyzer(const edm::ParameterSet& iConfig)
0019 : dqmStore_(edm::Service<DQMStore>().operator->()),
0020 fullconf_(iConfig),
0021 conf_(fullconf_.getParameter<edm::ParameterSet>("dEdxParameters")),
0022 doAllPlots_(conf_.getParameter<bool>("doAllPlots")),
0023 doDeDxPlots_(conf_.getParameter<bool>("doDeDxPlots")),
0024 genTriggerEventFlag_(new GenericTriggerEventFlag(
0025 conf_.getParameter<edm::ParameterSet>("genericTriggerEventPSet"), consumesCollector(), *this)) {
0026 trackInputTag_ = edm::InputTag(conf_.getParameter<std::string>("TracksForDeDx"));
0027 trackToken_ = consumes<reco::TrackCollection>(trackInputTag_);
0028
0029 dEdxInputList_ = conf_.getParameter<std::vector<std::string> >("deDxProducers");
0030 for (auto const& tag : dEdxInputList_) {
0031 dEdxTokenList_.push_back(consumes<reco::DeDxDataValueMap>(edm::InputTag(tag)));
0032 }
0033 }
0034
0035 dEdxAnalyzer::~dEdxAnalyzer() {
0036 if (genTriggerEventFlag_)
0037 delete genTriggerEventFlag_;
0038 }
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051 void dEdxAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) {
0052
0053 if (genTriggerEventFlag_->on())
0054 genTriggerEventFlag_->initRun(iRun, iSetup);
0055
0056
0057 std::string MEFolderName = conf_.getParameter<std::string>("FolderName");
0058
0059
0060 TrackHitMin = conf_.getParameter<double>("TrackHitMin");
0061 HIPdEdxMin = conf_.getParameter<double>("HIPdEdxMin");
0062 HighPtThreshold = conf_.getParameter<double>("HighPtThreshold");
0063
0064 dEdxK = conf_.getParameter<double>("dEdxK");
0065 dEdxC = conf_.getParameter<double>("dEdxC");
0066
0067 int dEdxNHitBin = conf_.getParameter<int>("dEdxNHitBin");
0068 double dEdxNHitMin = conf_.getParameter<double>("dEdxNHitMin");
0069 double dEdxNHitMax = conf_.getParameter<double>("dEdxNHitMax");
0070
0071 int dEdxBin = conf_.getParameter<int>("dEdxBin");
0072 double dEdxMin = conf_.getParameter<double>("dEdxMin");
0073 double dEdxMax = conf_.getParameter<double>("dEdxMax");
0074
0075 int dEdxHIPmassBin = conf_.getParameter<int>("dEdxHIPmassBin");
0076 double dEdxHIPmassMin = conf_.getParameter<double>("dEdxHIPmassMin");
0077 double dEdxHIPmassMax = conf_.getParameter<double>("dEdxHIPmassMax");
0078
0079 int dEdxMIPmassBin = conf_.getParameter<int>("dEdxMIPmassBin");
0080 double dEdxMIPmassMin = conf_.getParameter<double>("dEdxMIPmassMin");
0081 double dEdxMIPmassMax = conf_.getParameter<double>("dEdxMIPmassMax");
0082
0083 ibooker.setCurrentFolder(MEFolderName);
0084
0085
0086
0087
0088 if (doDeDxPlots_ || doAllPlots_) {
0089 for (unsigned int i = 0; i < dEdxInputList_.size(); i++) {
0090 ibooker.setCurrentFolder(MEFolderName + "/" + dEdxInputList_[i]);
0091 dEdxMEsVector.push_back(dEdxMEs());
0092
0093 histname = "MIP_dEdxPerTrack_";
0094 dEdxMEsVector[i].ME_MipDeDx = ibooker.book1D(histname, histname, dEdxBin, dEdxMin, dEdxMax);
0095 dEdxMEsVector[i].ME_MipDeDx->setAxisTitle("dEdx of each MIP Track (MeV/cm)");
0096 dEdxMEsVector[i].ME_MipDeDx->setAxisTitle("Number of Tracks", 2);
0097
0098 histname = "MIP_NumberOfdEdxHitsPerTrack_";
0099 dEdxMEsVector[i].ME_MipDeDxNHits = ibooker.book1D(histname, histname, dEdxNHitBin, dEdxNHitMin, dEdxNHitMax);
0100 dEdxMEsVector[i].ME_MipDeDxNHits->setAxisTitle("Number of dEdxHits of each MIP Track");
0101 dEdxMEsVector[i].ME_MipDeDxNHits->setAxisTitle("Number of Tracks", 2);
0102
0103 histname = "MIP_FractionOfSaturateddEdxHitsPerTrack_";
0104 dEdxMEsVector[i].ME_MipDeDxNSatHits = ibooker.book1D(histname, histname, 2 * dEdxNHitBin, 0, 1);
0105 dEdxMEsVector[i].ME_MipDeDxNSatHits->setAxisTitle("Fraction of Saturated dEdxHits of each MIP Track");
0106 dEdxMEsVector[i].ME_MipDeDxNSatHits->setAxisTitle("Number of Tracks", 2);
0107
0108 histname = "MIP_MassPerTrack_";
0109 dEdxMEsVector[i].ME_MipDeDxMass =
0110 ibooker.book1D(histname, histname, dEdxMIPmassBin, dEdxMIPmassMin, dEdxMIPmassMax);
0111 dEdxMEsVector[i].ME_MipDeDxMass->setAxisTitle("dEdx Mass of each MIP Track (GeV/c^{2})");
0112 dEdxMEsVector[i].ME_MipDeDxMass->setAxisTitle("Number of Tracks", 2);
0113
0114 histname = "HIP_MassPerTrack_";
0115 dEdxMEsVector[i].ME_HipDeDxMass =
0116 ibooker.book1D(histname, histname, dEdxHIPmassBin, dEdxHIPmassMin, dEdxHIPmassMax);
0117 dEdxMEsVector[i].ME_HipDeDxMass->setAxisTitle("dEdx Mass of each HIP Track (GeV/c^{2})");
0118 dEdxMEsVector[i].ME_HipDeDxMass->setAxisTitle("Number of Tracks", 2);
0119
0120 histname = "MIPOfHighPt_dEdxPerTrack_";
0121 dEdxMEsVector[i].ME_MipHighPtDeDx = ibooker.book1D(histname, histname, dEdxBin, dEdxMin, dEdxMax);
0122 dEdxMEsVector[i].ME_MipHighPtDeDx->setAxisTitle("dEdx of each MIP (of High pT) Track (MeV/cm)");
0123 dEdxMEsVector[i].ME_MipHighPtDeDx->setAxisTitle("Number of Tracks", 2);
0124
0125 histname = "MIPOfHighPt_NumberOfdEdxHitsPerTrack_";
0126 dEdxMEsVector[i].ME_MipHighPtDeDxNHits =
0127 ibooker.book1D(histname, histname, dEdxNHitBin, dEdxNHitMin, dEdxNHitMax);
0128 dEdxMEsVector[i].ME_MipHighPtDeDxNHits->setAxisTitle("Number of dEdxHits of each MIP (of High pT) Track");
0129 dEdxMEsVector[i].ME_MipHighPtDeDxNHits->setAxisTitle("Number of Tracks", 2);
0130 }
0131 }
0132 }
0133
0134 double dEdxAnalyzer::mass(double P, double I) {
0135 if (I - dEdxC < 0)
0136 return -1;
0137 return sqrt((I - dEdxC) / dEdxK) * P;
0138 }
0139
0140
0141
0142 void dEdxAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0143
0144 if (genTriggerEventFlag_->on() && !genTriggerEventFlag_->accept(iEvent, iSetup))
0145 return;
0146
0147 if (doDeDxPlots_ || doAllPlots_) {
0148 edm::Handle<reco::TrackCollection> trackCollectionHandle = iEvent.getHandle(trackToken_);
0149 if (!trackCollectionHandle.isValid())
0150 return;
0151
0152 for (unsigned int i = 0; i < dEdxInputList_.size(); i++) {
0153 edm::Handle<reco::DeDxDataValueMap> dEdxObjectHandle = iEvent.getHandle(dEdxTokenList_[i]);
0154 if (!dEdxObjectHandle.isValid())
0155 continue;
0156 const edm::ValueMap<reco::DeDxData> dEdxColl = *dEdxObjectHandle.product();
0157
0158 for (unsigned int t = 0; t < trackCollectionHandle->size(); t++) {
0159 reco::TrackRef track = reco::TrackRef(trackCollectionHandle, t);
0160
0161 if (track->quality(reco::TrackBase::highPurity)) {
0162
0163 if (track->pt() >= 5.0 && track->numberOfValidHits() > TrackHitMin) {
0164 dEdxMEsVector[i].ME_MipDeDx->Fill(dEdxColl[track].dEdx());
0165 dEdxMEsVector[i].ME_MipDeDxNHits->Fill(dEdxColl[track].numberOfMeasurements());
0166 if (dEdxColl[track].numberOfMeasurements() != 0)
0167 dEdxMEsVector[i].ME_MipDeDxNSatHits->Fill((1.0 * dEdxColl[track].numberOfSaturatedMeasurements()) /
0168 dEdxColl[track].numberOfMeasurements());
0169 dEdxMEsVector[i].ME_MipDeDxMass->Fill(mass(track->p(), dEdxColl[track].dEdx()));
0170
0171 if (track->pt() >= HighPtThreshold) {
0172 dEdxMEsVector[i].ME_MipHighPtDeDx->Fill(dEdxColl[track].dEdx());
0173 dEdxMEsVector[i].ME_MipHighPtDeDxNHits->Fill(dEdxColl[track].numberOfMeasurements());
0174 }
0175
0176
0177 } else if (track->pt() < 2 && dEdxColl[track].dEdx() > HIPdEdxMin) {
0178 dEdxMEsVector[i].ME_HipDeDxMass->Fill(mass(track->p(), dEdxColl[track].dEdx()));
0179 }
0180 }
0181 }
0182 }
0183 }
0184 }
0185
0186
0187 void dEdxAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0188
0189
0190 edm::ParameterSetDescription desc;
0191 desc.setUnknown();
0192 descriptions.addDefault(desc);
0193 }
0194
0195 DEFINE_FWK_MODULE(dEdxAnalyzer);