TrackingMonitor/src/dEdxAnalyzer.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author Loic Quertenmont
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 // -- BeginRun
0042 //---------------------------------------------------------------------------------//
0043 void dEdxAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
0044 {
0045
0046   // Initialize the GenericTriggerEventFlag
0047   if ( genTriggerEventFlag_->on() ) genTriggerEventFlag_->initRun( iRun, iSetup );
0048 }
0049 */
0050
0051 void dEdxAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) {
0052   // Initialize the GenericTriggerEventFlag
0053   if (genTriggerEventFlag_->on())
0054     genTriggerEventFlag_->initRun(iRun, iSetup);
0055
0056   // parameters from the configuration
0057   std::string MEFolderName = conf_.getParameter<std::string>("FolderName");
0058
0059   // get binning from the configuration
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   // book the Hit Property histograms
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 // -- Analyse
0141 // ---------------------------------------------------------------------------------//
0142 void dEdxAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0143   // Filter out events if Trigger Filtering is requested
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           //MIPs
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             //HighlyIonizing particles
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 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0187 void dEdxAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0188   //The following says we do not know what parameters are allowed so do no validation
0189   // Please change this to state exactly what you do use, even if it is no parameters
0190   edm::ParameterSetDescription desc;
0191   desc.setUnknown();
0192   descriptions.addDefault(desc);
0193 }
0194
0195 DEFINE_FWK_MODULE(dEdxAnalyzer);