Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-25 02:44:10 

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 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 // Constructor
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   //now do what ever initialization is needed
0113   theTrackCollectionToken = consumes<reco::TrackCollection>(TrackCollectionTag_);
0114   theVertexCollectionToken = consumes<reco::VertexCollection>(PVtxCollectionTag_);
0115 }
0116 
0117 //
0118 // member functions
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   // Initialize Root-tuple variables
0131   //=======================================================
0132 
0133   SetVarToZero();
0134 
0135   //=======================================================
0136   // Retrieve the Track information
0137   //=======================================================
0138 
0139   // Declare the handle for the vertex collection
0140   const auto& vertexHandle = iEvent.getHandle(theVertexCollectionToken);
0141 
0142   // Check if the handle is valid
0143   if (!vertexHandle.isValid()) {
0144     edm::LogError("LhcTrackAnalyzer") << "Vertex collection not found or invalid.";
0145     return;  // Early return if the vertex collection is not valid
0146   }
0147 
0148   // Retrieve the actual product from the handle
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   // Check if the handle is valid
0173   if (!tracksHandle.isValid()) {
0174     edm::LogError("LhcTrackAnalyzer") << "Tracks collection not found or invalid.";
0175     return;  // Early return if the vertex collection is not valid
0176   }
0177 
0178   // Retrieve the actual product from the handle
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     // Loop over vertices to find the closest one in dz
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       // Compute dxy and dz with respect to the closest vertex
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   }  //end loop on tracks
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   // Define TTree for output
0327   rootFile_ = new TFile(filename_.c_str(), "recreate");
0328   rootTree_ = new TTree("tree", "Lhc Track tree");
0329 
0330   // Track Paramters
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 //define this as a plug-in
0416 DEFINE_FWK_MODULE(LhcTrackAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!