File indexing completed on 2024-06-25 02:44:10
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 dzPV_[nMaxtracks_];
0086 double dxyPV_[nMaxtracks_];
0087 double dzSig_[nMaxtracks_];
0088 double dxySig_[nMaxtracks_];
0089 double dzPVSig_[nMaxtracks_];
0090 double dxyPVSig_[nMaxtracks_];
0091 double dzError_[nMaxtracks_];
0092 double dxyError_[nMaxtracks_];
0093 double xPCA_[nMaxtracks_];
0094 double yPCA_[nMaxtracks_];
0095 double zPCA_[nMaxtracks_];
0096 int trkAlgo_[nMaxtracks_];
0097 int trkQuality_[nMaxtracks_];
0098 int isHighPurity_[nMaxtracks_];
0099 int validhits_[nMaxtracks_][7];
0100 bool goodbx_;
0101 bool goodvtx_;
0102 };
0103
0104
0105
0106 LhcTrackAnalyzer::LhcTrackAnalyzer(const edm::ParameterSet& iConfig)
0107 : TrackCollectionTag_(iConfig.getParameter<edm::InputTag>("TrackCollectionTag")),
0108 PVtxCollectionTag_(iConfig.getParameter<edm::InputTag>("PVtxCollectionTag")),
0109 debug_(iConfig.getParameter<bool>("Debug")),
0110 acceptedBX_(iConfig.getParameter<std::vector<unsigned int>>("acceptedBX")),
0111 filename_(iConfig.getParameter<std::string>("OutputFileName")) {
0112
0113 theTrackCollectionToken = consumes<reco::TrackCollection>(TrackCollectionTag_);
0114 theVertexCollectionToken = consumes<reco::VertexCollection>(PVtxCollectionTag_);
0115 }
0116
0117
0118
0119
0120
0121
0122 void LhcTrackAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
0123
0124 {
0125 using namespace edm;
0126 using namespace reco;
0127 using namespace std;
0128
0129
0130
0131
0132
0133 SetVarToZero();
0134
0135
0136
0137
0138
0139
0140 const auto& vertexHandle = iEvent.getHandle(theVertexCollectionToken);
0141
0142
0143 if (!vertexHandle.isValid()) {
0144 edm::LogError("LhcTrackAnalyzer") << "Vertex collection not found or invalid.";
0145 return;
0146 }
0147
0148
0149 const reco::VertexCollection& vertices = *vertexHandle;
0150
0151 const auto& vtx = vertices.front();
0152 if (vtx.isFake()) {
0153 goodvtx_ = false;
0154 } else {
0155 goodvtx_ = true;
0156 }
0157
0158 int bx = iEvent.bunchCrossing();
0159 if (acceptedBX_.empty()) {
0160 goodbx_ = true;
0161 } else {
0162 if (std::find(acceptedBX_.begin(), acceptedBX_.end(), bx) != acceptedBX_.end()) {
0163 goodbx_ = true;
0164 }
0165 }
0166
0167 run_ = iEvent.id().run();
0168 event_ = iEvent.id().event();
0169
0170 const auto& tracksHandle = iEvent.getHandle(theTrackCollectionToken);
0171
0172
0173 if (!tracksHandle.isValid()) {
0174 edm::LogError("LhcTrackAnalyzer") << "Tracks collection not found or invalid.";
0175 return;
0176 }
0177
0178
0179 const reco::TrackCollection& tracks = *tracksHandle;
0180
0181 if (debug_) {
0182 edm::LogInfo("LhcTrackAnalyzer") << "LhcTrackAnalyzer::analyze() looping over " << tracks.size() << "tracks."
0183 << endl;
0184 }
0185
0186 for (const auto& track : tracks) {
0187 if (nTracks_ >= nMaxtracks_) {
0188 edm::LogWarning("LhcTrackAnalyzer")
0189 << " LhcTrackAnalyzer::analyze() : Warning - Run " << run_ << " Event " << event_
0190 << "\tNumber of tracks: " << tracks.size() << " , greater than " << nMaxtracks_ << std::endl;
0191 continue;
0192 }
0193 pt_[nTracks_] = track.pt();
0194 eta_[nTracks_] = track.eta();
0195 phi_[nTracks_] = track.phi();
0196 chi2_[nTracks_] = track.chi2();
0197 chi2ndof_[nTracks_] = track.normalizedChi2();
0198 charge_[nTracks_] = track.charge();
0199 qoverp_[nTracks_] = track.qoverp();
0200 dz_[nTracks_] = track.dz();
0201 dxy_[nTracks_] = track.dxy();
0202 dzError_[nTracks_] = track.dzError();
0203 dxyError_[nTracks_] = track.dxyError();
0204 dzSig_[nTracks_] = track.dz() / track.dzError();
0205 dxySig_[nTracks_] = track.dxy() / track.dxyError();
0206
0207 const reco::Vertex* closestVertex = nullptr;
0208 float minDz = std::numeric_limits<float>::max();
0209
0210
0211 for (const auto& vertex : vertices) {
0212 float dz = fabs(track.dz(vertex.position()));
0213 if (dz < minDz) {
0214 minDz = dz;
0215 closestVertex = &vertex;
0216 }
0217 }
0218
0219 float dzPV{-999.};
0220 float dxyPV{-999.};
0221 if (closestVertex) {
0222
0223 dzPV = track.dz(closestVertex->position());
0224 dxyPV = track.dxy(closestVertex->position());
0225 }
0226
0227 dzPV_[nTracks_] = dzPV;
0228 dxyPV_[nTracks_] = dxyPV;
0229
0230 dzPVSig_[nTracks_] = dzPV / track.dzError();
0231 dxyPVSig_[nTracks_] = dxyPV / track.dxyError();
0232
0233 xPCA_[nTracks_] = track.vertex().x();
0234 yPCA_[nTracks_] = track.vertex().y();
0235 zPCA_[nTracks_] = track.vertex().z();
0236 validhits_[nTracks_][0] = track.numberOfValidHits();
0237 validhits_[nTracks_][1] = track.hitPattern().numberOfValidPixelBarrelHits();
0238 validhits_[nTracks_][2] = track.hitPattern().numberOfValidPixelEndcapHits();
0239 validhits_[nTracks_][3] = track.hitPattern().numberOfValidStripTIBHits();
0240 validhits_[nTracks_][4] = track.hitPattern().numberOfValidStripTIDHits();
0241 validhits_[nTracks_][5] = track.hitPattern().numberOfValidStripTOBHits();
0242 validhits_[nTracks_][6] = track.hitPattern().numberOfValidStripTECHits();
0243
0244 int myalgo = -88;
0245 if (track.algo() == reco::TrackBase::undefAlgorithm) {
0246 myalgo = 0;
0247 } else if (track.algo() == reco::TrackBase::ctf) {
0248 myalgo = 1;
0249 } else if (track.algo() == reco::TrackBase::duplicateMerge) {
0250 myalgo = 2;
0251 } else if (track.algo() == reco::TrackBase::cosmics) {
0252 myalgo = 3;
0253 } else if (track.algo() == reco::TrackBase::initialStep) {
0254 myalgo = 4;
0255 } else if (track.algo() == reco::TrackBase::lowPtTripletStep) {
0256 myalgo = 5;
0257 } else if (track.algo() == reco::TrackBase::pixelPairStep) {
0258 myalgo = 6;
0259 } else if (track.algo() == reco::TrackBase::detachedTripletStep) {
0260 myalgo = 7;
0261 } else if (track.algo() == reco::TrackBase::mixedTripletStep) {
0262 myalgo = 8;
0263 } else if (track.algo() == reco::TrackBase::pixelLessStep) {
0264 myalgo = 9;
0265 } else if (track.algo() == reco::TrackBase::tobTecStep) {
0266 myalgo = 10;
0267 } else if (track.algo() == reco::TrackBase::jetCoreRegionalStep) {
0268 myalgo = 11;
0269 } else if (track.algo() == reco::TrackBase::muonSeededStepInOut) {
0270 myalgo = 13;
0271 } else if (track.algo() == reco::TrackBase::muonSeededStepOutIn) {
0272 myalgo = 14;
0273 } else if (track.algo() == reco::TrackBase::highPtTripletStep) {
0274 myalgo = 22;
0275 } else if (track.algo() == reco::TrackBase::lowPtQuadStep) {
0276 myalgo = 23;
0277 } else if (track.algo() == reco::TrackBase::detachedQuadStep) {
0278 myalgo = 24;
0279 } else if (track.algo() == reco::TrackBase::displacedGeneralStep) {
0280 myalgo = 25;
0281 } else if (track.algo() == reco::TrackBase::displacedRegionalStep) {
0282 myalgo = 26;
0283 } else if (track.algo() == reco::TrackBase::hltIter0) {
0284 myalgo = 31;
0285 } else {
0286 myalgo = reco::TrackBase::algoSize;
0287 edm::LogWarning("LhcTrackAnalyzer")
0288 << "LhcTrackAnalyzer does not support all types of tracks, encountered one from algo "
0289 << reco::TrackBase::algoName(track.algo());
0290 }
0291 trkAlgo_[nTracks_] = myalgo;
0292
0293 int myquality = -99;
0294 if (track.quality(reco::TrackBase::undefQuality))
0295 myquality = -1;
0296 if (track.quality(reco::TrackBase::loose))
0297 myquality = 0;
0298 if (track.quality(reco::TrackBase::tight))
0299 myquality = 1;
0300 if (track.quality(reco::TrackBase::highPurity))
0301 myquality = 2;
0302 trkQuality_[nTracks_] = myquality;
0303
0304 if (track.quality(reco::TrackBase::highPurity))
0305 isHighPurity_[nTracks_] = 1;
0306 else
0307 isHighPurity_[nTracks_] = 0;
0308 nTracks_++;
0309
0310 }
0311
0312 for (int d = 0; d < nTracks_; ++d) {
0313 if (abs(trkQuality_[d]) > 5)
0314 edm::LogInfo("LhcTrackAnalyzer") << "MYQUALITY!!! " << trkQuality_[d] << " at track # " << d << "/" << nTracks_
0315 << endl;
0316 }
0317
0318 rootTree_->Fill();
0319 }
0320
0321
0322 void LhcTrackAnalyzer::beginJob()
0323
0324 {
0325 edm::LogInfo("beginJob") << "Begin Job" << std::endl;
0326
0327 rootFile_ = new TFile(filename_.c_str(), "recreate");
0328 rootTree_ = new TTree("tree", "Lhc Track tree");
0329
0330
0331 rootTree_->Branch("run", &run_, "run/I");
0332 rootTree_->Branch("event", &event_, "event/I");
0333 rootTree_->Branch("goodbx", &goodbx_, "goodbx/O");
0334 rootTree_->Branch("goodvtx", &goodvtx_, "goodvtx/O");
0335 rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
0336 rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
0337 rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
0338 rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
0339 rootTree_->Branch("chi2", &chi2_, "chi2[nTracks]/D");
0340 rootTree_->Branch("chi2ndof", &chi2ndof_, "chi2ndof[nTracks]/D");
0341 rootTree_->Branch("charge", &charge_, "charge[nTracks]/I");
0342 rootTree_->Branch("qoverp", &qoverp_, "qoverp[nTracks]/D");
0343 rootTree_->Branch("dz", &dz_, "dz[nTracks]/D");
0344 rootTree_->Branch("dxy", &dxy_, "dxy[nTracks]/D");
0345 rootTree_->Branch("dzPV", &dzPV_, "dzPV[nTracks]/D");
0346 rootTree_->Branch("dxyPV", &dxy_, "dxyPV[nTracks]/D");
0347 rootTree_->Branch("dzError", &dzError_, "dzError[nTracks]/D");
0348 rootTree_->Branch("dxyError", &dxyError_, "dxyError[nTracks]/D");
0349 rootTree_->Branch("dzSig", &dzSig_, "dzSig[nTracks]/D");
0350 rootTree_->Branch("dxySig", &dxySig_, "dxySig[nTracks]/D");
0351 rootTree_->Branch("dzPVSig", &dzPVSig_, "dzPVSig[nTracks]/D");
0352 rootTree_->Branch("dxyPVSig", &dxyPVSig_, "dxyPVSig[nTracks]/D");
0353 rootTree_->Branch("xPCA", &xPCA_, "xPCA[nTracks]/D");
0354 rootTree_->Branch("yPCA", &yPCA_, "yPCA[nTracks]/D");
0355 rootTree_->Branch("zPCA", &zPCA_, "zPCA[nTracks]/D");
0356 rootTree_->Branch("isHighPurity", &isHighPurity_, "isHighPurity[nTracks]/I");
0357 rootTree_->Branch("trkQuality", &trkQuality_, "trkQuality[nTracks]/I");
0358 rootTree_->Branch("trkAlgo", &trkAlgo_, "trkAlgo[nTracks]/I");
0359 rootTree_->Branch("nValidHits", &validhits_, "nValidHits[nTracks][7]/I");
0360 }
0361
0362
0363 void LhcTrackAnalyzer::endJob()
0364
0365 {
0366 if (rootFile_) {
0367 rootFile_->Write();
0368 rootFile_->Close();
0369 }
0370 }
0371
0372
0373 void LhcTrackAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
0374
0375 {
0376 edm::ParameterSetDescription desc;
0377 desc.setComment("Ntuplizer for LHC tracks");
0378 desc.add<edm::InputTag>("TrackCollectionTag", edm::InputTag("ALCARECOTkAlMinBias"));
0379 desc.add<edm::InputTag>("PVtxCollectionTag", edm::InputTag("offlinePrimaryVertices"));
0380 desc.add<bool>("Debug", false);
0381 desc.add<std::vector<unsigned int>>("acceptedBX", {});
0382 desc.add<std::string>("OutputFileName", "LhcTrackAnalyzer_Output_default.root");
0383 descriptions.addWithDefaultLabel(desc);
0384 }
0385
0386
0387 void LhcTrackAnalyzer::SetVarToZero()
0388
0389 {
0390 run_ = -1;
0391 event_ = -99;
0392 nTracks_ = 0;
0393 for (int i = 0; i < nMaxtracks_; ++i) {
0394 pt_[i] = 0;
0395 eta_[i] = 0;
0396 phi_[i] = 0;
0397 chi2_[i] = 0;
0398 chi2ndof_[i] = 0;
0399 charge_[i] = 0;
0400 qoverp_[i] = 0;
0401 dz_[i] = 0;
0402 dxy_[i] = 0;
0403 xPCA_[i] = 0;
0404 yPCA_[i] = 0;
0405 zPCA_[i] = 0;
0406 trkQuality_[i] = 0;
0407 trkAlgo_[i] = -1;
0408 isHighPurity_[i] = -3;
0409 for (int j = 0; j < 7; j++) {
0410 validhits_[nTracks_][j] = -1 * j;
0411 }
0412 }
0413 }
0414
0415
0416 DEFINE_FWK_MODULE(LhcTrackAnalyzer);