File indexing completed on 2021-08-17 23:10:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include <memory>
0017 #include <string>
0018 #include <vector>
0019
0020
0021 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0022 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0025 #include "DataFormats/VertexReco/interface/Vertex.h"
0026 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0027 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/EventSetup.h"
0030 #include "FWCore/Framework/interface/Frameworkfwd.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "FWCore/ServiceRegistry/interface/Service.h"
0034
0035
0036 #include "TFile.h"
0037 #include "TTree.h"
0038
0039
0040
0041
0042
0043 class LhcTrackAnalyzer : public edm::one::EDAnalyzer<> {
0044 public:
0045 explicit LhcTrackAnalyzer(const edm::ParameterSet&);
0046 ~LhcTrackAnalyzer() override = default;
0047 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0048
0049 private:
0050 void beginJob() override;
0051 void analyze(const edm::Event&, const edm::EventSetup&) override;
0052 void endJob() override;
0053
0054
0055 edm::InputTag TrackCollectionTag_;
0056 edm::InputTag PVtxCollectionTag_;
0057 bool debug_;
0058 std::vector<unsigned int> acceptedBX_;
0059
0060 edm::EDGetTokenT<reco::TrackCollection> theTrackCollectionToken;
0061 edm::EDGetTokenT<reco::VertexCollection> theVertexCollectionToken;
0062
0063
0064 std::string filename_;
0065 TFile* rootFile_;
0066 TTree* rootTree_;
0067
0068
0069
0070 void SetVarToZero();
0071
0072 static constexpr int nMaxtracks_ = 3000;
0073 int nTracks_;
0074 int run_;
0075 int event_;
0076 double pt_[nMaxtracks_];
0077 double eta_[nMaxtracks_];
0078 double phi_[nMaxtracks_];
0079 double chi2_[nMaxtracks_];
0080 double chi2ndof_[nMaxtracks_];
0081 int charge_[nMaxtracks_];
0082 double qoverp_[nMaxtracks_];
0083 double dz_[nMaxtracks_];
0084 double dxy_[nMaxtracks_];
0085 double xPCA_[nMaxtracks_];
0086 double yPCA_[nMaxtracks_];
0087 double zPCA_[nMaxtracks_];
0088 int trkAlgo_[nMaxtracks_];
0089 int trkQuality_[nMaxtracks_];
0090 int isHighPurity_[nMaxtracks_];
0091 int validhits_[nMaxtracks_][7];
0092 bool goodbx_;
0093 bool goodvtx_;
0094 };
0095
0096
0097
0098 LhcTrackAnalyzer::LhcTrackAnalyzer(const edm::ParameterSet& iConfig)
0099 : TrackCollectionTag_(iConfig.getParameter<edm::InputTag>("TrackCollectionTag")),
0100 PVtxCollectionTag_(iConfig.getParameter<edm::InputTag>("PVtxCollectionTag")),
0101 debug_(iConfig.getParameter<bool>("Debug")),
0102 acceptedBX_(iConfig.getParameter<std::vector<unsigned int>>("acceptedBX")),
0103 filename_(iConfig.getParameter<std::string>("OutputFileName")) {
0104
0105 theTrackCollectionToken = consumes<reco::TrackCollection>(TrackCollectionTag_);
0106 theVertexCollectionToken = consumes<reco::VertexCollection>(PVtxCollectionTag_);
0107 }
0108
0109
0110
0111
0112
0113
0114 void LhcTrackAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
0115
0116 {
0117 using namespace edm;
0118 using namespace reco;
0119 using namespace std;
0120
0121
0122
0123
0124
0125 SetVarToZero();
0126
0127
0128
0129
0130
0131 const auto& vertices = iEvent.get(theVertexCollectionToken);
0132 const auto& vtx = vertices.front();
0133 if (vtx.isFake()) {
0134 goodvtx_ = false;
0135 } else {
0136 goodvtx_ = true;
0137 }
0138
0139 int bx = iEvent.bunchCrossing();
0140 if (acceptedBX_.empty()) {
0141 goodbx_ = true;
0142 } else {
0143 if (std::find(acceptedBX_.begin(), acceptedBX_.end(), bx) != acceptedBX_.end()) {
0144 goodbx_ = true;
0145 }
0146 }
0147
0148 run_ = iEvent.id().run();
0149 event_ = iEvent.id().event();
0150
0151 const auto& tracks = iEvent.get(theTrackCollectionToken);
0152 if (debug_) {
0153 edm::LogInfo("LhcTrackAnalyzer") << "LhcTrackAnalyzer::analyze() looping over " << tracks.size() << "tracks."
0154 << endl;
0155 }
0156
0157 for (const auto& track : tracks) {
0158 if (nTracks_ >= nMaxtracks_) {
0159 edm::LogWarning("LhcTrackAnalyzer")
0160 << " LhcTrackAnalyzer::analyze() : Warning - Run " << run_ << " Event " << event_
0161 << "\tNumber of tracks: " << tracks.size() << " , greater than " << nMaxtracks_ << std::endl;
0162 continue;
0163 }
0164 pt_[nTracks_] = track.pt();
0165 eta_[nTracks_] = track.eta();
0166 phi_[nTracks_] = track.phi();
0167 chi2_[nTracks_] = track.chi2();
0168 chi2ndof_[nTracks_] = track.normalizedChi2();
0169 charge_[nTracks_] = track.charge();
0170 qoverp_[nTracks_] = track.qoverp();
0171 dz_[nTracks_] = track.dz();
0172 dxy_[nTracks_] = track.dxy();
0173 xPCA_[nTracks_] = track.vertex().x();
0174 yPCA_[nTracks_] = track.vertex().y();
0175 zPCA_[nTracks_] = track.vertex().z();
0176 validhits_[nTracks_][0] = track.numberOfValidHits();
0177 validhits_[nTracks_][1] = track.hitPattern().numberOfValidPixelBarrelHits();
0178 validhits_[nTracks_][2] = track.hitPattern().numberOfValidPixelEndcapHits();
0179 validhits_[nTracks_][3] = track.hitPattern().numberOfValidStripTIBHits();
0180 validhits_[nTracks_][4] = track.hitPattern().numberOfValidStripTIDHits();
0181 validhits_[nTracks_][5] = track.hitPattern().numberOfValidStripTOBHits();
0182 validhits_[nTracks_][6] = track.hitPattern().numberOfValidStripTECHits();
0183
0184 int myalgo = -88;
0185 if (track.algo() == reco::TrackBase::undefAlgorithm) {
0186 myalgo = 0;
0187 } else if (track.algo() == reco::TrackBase::ctf) {
0188 myalgo = 1;
0189 } else if (track.algo() == reco::TrackBase::duplicateMerge) {
0190 myalgo = 2;
0191 } else if (track.algo() == reco::TrackBase::cosmics) {
0192 myalgo = 3;
0193 } else if (track.algo() == reco::TrackBase::initialStep) {
0194 myalgo = 4;
0195 } else if (track.algo() == reco::TrackBase::lowPtTripletStep) {
0196 myalgo = 5;
0197 } else if (track.algo() == reco::TrackBase::pixelPairStep) {
0198 myalgo = 6;
0199 } else if (track.algo() == reco::TrackBase::detachedTripletStep) {
0200 myalgo = 7;
0201 } else if (track.algo() == reco::TrackBase::mixedTripletStep) {
0202 myalgo = 8;
0203 } else if (track.algo() == reco::TrackBase::pixelLessStep) {
0204 myalgo = 9;
0205 } else if (track.algo() == reco::TrackBase::tobTecStep) {
0206 myalgo = 10;
0207 } else if (track.algo() == reco::TrackBase::jetCoreRegionalStep) {
0208 myalgo = 11;
0209 } else if (track.algo() == reco::TrackBase::muonSeededStepInOut) {
0210 myalgo = 13;
0211 } else if (track.algo() == reco::TrackBase::muonSeededStepOutIn) {
0212 myalgo = 14;
0213 } else if (track.algo() == reco::TrackBase::highPtTripletStep) {
0214 myalgo = 22;
0215 } else if (track.algo() == reco::TrackBase::lowPtQuadStep) {
0216 myalgo = 23;
0217 } else if (track.algo() == reco::TrackBase::detachedQuadStep) {
0218 myalgo = 24;
0219 } else {
0220 myalgo = 25;
0221 edm::LogWarning("LhcTrackAnalyzer")
0222 << "LhcTrackAnalyzer does not support all types of tracks, encountered one from algo "
0223 << reco::TrackBase::algoName(track.algo());
0224 }
0225 trkAlgo_[nTracks_] = myalgo;
0226
0227 int myquality = -99;
0228 if (track.quality(reco::TrackBase::undefQuality))
0229 myquality = -1;
0230 if (track.quality(reco::TrackBase::loose))
0231 myquality = 0;
0232 if (track.quality(reco::TrackBase::tight))
0233 myquality = 1;
0234 if (track.quality(reco::TrackBase::highPurity))
0235 myquality = 2;
0236 trkQuality_[nTracks_] = myquality;
0237
0238 if (track.quality(reco::TrackBase::highPurity))
0239 isHighPurity_[nTracks_] = 1;
0240 else
0241 isHighPurity_[nTracks_] = 0;
0242 nTracks_++;
0243
0244 }
0245
0246 for (int d = 0; d < nTracks_; ++d) {
0247 if (abs(trkQuality_[d]) > 5)
0248 edm::LogInfo("LhcTrackAnalyzer") << "MYQUALITY!!! " << trkQuality_[d] << " at track # " << d << "/" << nTracks_
0249 << endl;
0250 }
0251
0252 rootTree_->Fill();
0253 }
0254
0255
0256 void LhcTrackAnalyzer::beginJob()
0257
0258 {
0259 edm::LogInfo("beginJob") << "Begin Job" << std::endl;
0260
0261 rootFile_ = new TFile(filename_.c_str(), "recreate");
0262 rootTree_ = new TTree("tree", "Lhc Track tree");
0263
0264
0265 rootTree_->Branch("run", &run_, "run/I");
0266 rootTree_->Branch("event", &event_, "event/I");
0267 rootTree_->Branch("goodbx", &goodbx_, "goodbx/O");
0268 rootTree_->Branch("goodvtx", &goodvtx_, "goodvtx/O");
0269 rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
0270 rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
0271 rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
0272 rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
0273 rootTree_->Branch("chi2", &chi2_, "chi2[nTracks]/D");
0274 rootTree_->Branch("chi2ndof", &chi2ndof_, "chi2ndof[nTracks]/D");
0275 rootTree_->Branch("charge", &charge_, "charge[nTracks]/I");
0276 rootTree_->Branch("qoverp", &qoverp_, "qoverp[nTracks]/D");
0277 rootTree_->Branch("dz", &dz_, "dz[nTracks]/D");
0278 rootTree_->Branch("dxy", &dxy_, "dxy[nTracks]/D");
0279 rootTree_->Branch("xPCA", &xPCA_, "xPCA[nTracks]/D");
0280 rootTree_->Branch("yPCA", &yPCA_, "yPCA[nTracks]/D");
0281 rootTree_->Branch("zPCA", &zPCA_, "zPCA[nTracks]/D");
0282 rootTree_->Branch("isHighPurity", &isHighPurity_, "isHighPurity[nTracks]/I");
0283 rootTree_->Branch("trkQuality", &trkQuality_, "trkQuality[nTracks]/I");
0284 rootTree_->Branch("trkAlgo", &trkAlgo_, "trkAlgo[nTracks]/I");
0285 rootTree_->Branch("nValidHits", &validhits_, "nValidHits[nTracks][7]/I");
0286 }
0287
0288
0289 void LhcTrackAnalyzer::endJob()
0290
0291 {
0292 if (rootFile_) {
0293 rootFile_->Write();
0294 rootFile_->Close();
0295 }
0296 }
0297
0298
0299 void LhcTrackAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
0300
0301 {
0302 edm::ParameterSetDescription desc;
0303 desc.setComment("Ntuplizer for LHC tracks");
0304 desc.add<edm::InputTag>("TrackCollectionTag", edm::InputTag("ALCARECOTkAlMinBias"));
0305 desc.add<edm::InputTag>("PVtxCollectionTag", edm::InputTag("offlinePrimaryVertices"));
0306 desc.add<bool>("Debug", false);
0307 desc.add<std::vector<unsigned int>>("acceptedBX", {});
0308 desc.add<std::string>("OutputFileName", "LhcTrackAnalyzer_Output_default.root");
0309 descriptions.addWithDefaultLabel(desc);
0310 }
0311
0312
0313 void LhcTrackAnalyzer::SetVarToZero()
0314
0315 {
0316 run_ = -1;
0317 event_ = -99;
0318 nTracks_ = 0;
0319 for (int i = 0; i < nMaxtracks_; ++i) {
0320 pt_[i] = 0;
0321 eta_[i] = 0;
0322 phi_[i] = 0;
0323 chi2_[i] = 0;
0324 chi2ndof_[i] = 0;
0325 charge_[i] = 0;
0326 qoverp_[i] = 0;
0327 dz_[i] = 0;
0328 dxy_[i] = 0;
0329 xPCA_[i] = 0;
0330 yPCA_[i] = 0;
0331 zPCA_[i] = 0;
0332 trkQuality_[i] = 0;
0333 trkAlgo_[i] = -1;
0334 isHighPurity_[i] = -3;
0335 for (int j = 0; j < 7; j++) {
0336 validhits_[nTracks_][j] = -1 * j;
0337 }
0338 }
0339 }
0340
0341
0342 DEFINE_FWK_MODULE(LhcTrackAnalyzer);
