Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-17 23:10:19

0001 // -*- C++ -*-
0002 //
0003 // Package:    LhcTrackAnalyzer
0004 // Class:      LhcTrackAnalyzer
0005 //
0006 /**\class LhcTrackAnalyzer LhcTrackAnalyzer.cc Alignment/HIPAlignmentAlgorithm/plugins/LhcTrackAnalyzer.cc
0007 
0008    Originally written by M.Musich
0009    Expanded by A. Bonato
0010 
0011    Description: Ntuplizer for collision tracks
0012 */
0013 //
0014 
0015 // system include files
0016 #include <memory>
0017 #include <string>
0018 #include <vector>
0019 
0020 // user include files
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 // ROOT includes
0036 #include "TFile.h"
0037 #include "TTree.h"
0038 
0039 //
0040 // class decleration
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   // ----------member data ---------------------------
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   // Output
0064   std::string filename_;
0065   TFile* rootFile_;
0066   TTree* rootTree_;
0067 
0068   // Root-Tuple variables :
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 // Constructor
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   //now do what ever initialization is needed
0105   theTrackCollectionToken = consumes<reco::TrackCollection>(TrackCollectionTag_);
0106   theVertexCollectionToken = consumes<reco::VertexCollection>(PVtxCollectionTag_);
0107 }
0108 
0109 //
0110 // member functions
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   // Initialize Root-tuple variables
0123   //=======================================================
0124 
0125   SetVarToZero();
0126 
0127   //=======================================================
0128   // Retrieve the Track information
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   }  //end loop on tracks
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   // Define TTree for output
0261   rootFile_ = new TFile(filename_.c_str(), "recreate");
0262   rootTree_ = new TTree("tree", "Lhc Track tree");
0263 
0264   // Track Paramters
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 //define this as a plug-in
0342 DEFINE_FWK_MODULE(LhcTrackAnalyzer);