TrackingMonitor/src/dEdxHitAnalyzer.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author Loic Quertenmont
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   // parameters from the configuration
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 // -- BeginRun
0063 //---------------------------------------------------------------------------------//
0064 void dEdxHitAnalyzer::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0065   // Initialize the GenericTriggerEventFlag
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   // book the Hit Property histograms
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);  //save charge
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 // -- Analyse
0127 // ---------------------------------------------------------------------------------//
0128 void dEdxHitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0129   // Filter out events if Trigger Filtering is requested
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 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0172 void dEdxHitAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0173   //The following says we do not know what parameters are allowed so do no validation
0174   // Please change this to state exactly what you do use, even if it is no parameters
0175   edm::ParameterSetDescription desc;
0176   desc.setUnknown();
0177   descriptions.addDefault(desc);
0178 }
0179
0180 DEFINE_FWK_MODULE(dEdxHitAnalyzer);