Tracking/test/testGeneralTracks.cc

0001 // user include files
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006
0007 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0008 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0009 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0010 #include "DataFormats/DetId/interface/DetId.h"
0011 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0012
0013 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0014
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019
0020 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0021 #include "FastSimulation/Event/interface/FSimEvent.h"
0022 #include "FastSimulation/Event/interface/FSimTrack.h"
0023 #include "FastSimulation/Event/interface/FSimVertex.h"
0024 //#include "FastSimulation/Tracking/interface/TrajectorySeedHitCandidate.h"
0025
0026 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0027 #include "DQMServices/Core/interface/DQMStore.h"
0028 #include "FWCore/ServiceRegistry/interface/Service.h"
0029 #include <vector>
0030 #include <string>
0031 //#include "TFile.h"
0032 //#include "TTree.h"
0033 //#include "TProcessID.h"
0034
0035 class testGeneralTracks : public DQMEDAnalyzer {
0036 public:
0037   explicit testGeneralTracks(const edm::ParameterSet&);
0038   ~testGeneralTracks();
0039
0040   virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
0041   virtual void dqmBeginRun(edm::Run const&, edm::EventSetup const&) override;
0042   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0043
0044 private:
0045   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken;
0046   const edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord> pdtToken;
0047   const edm::EDGetTokenT<std::vector<SimTrack>> tokSimTrack_;
0048   const edm::EDGetTokenT<std::vector<SimVertex>> tokSimVertex_;
0049
0050   // See RecoParticleFlow/PFProducer/interface/PFProducer.h
0051   edm::ParameterSet particleFilter_;
0052   std::vector<edm::InputTag> allTracks_;
0053   std::vector<edm::EDGetTokenT<reco::TrackCollection>> allTrackToken_;
0054   bool saveNU;
0055   std::vector<FSimEvent*> mySimEvent;
0056   std::string simModuleLabel_;
0057   // Histograms
0058   std::vector<MonitorElement*> h0;
0059   MonitorElement* genTracksvsEtaP;
0060   std::vector<MonitorElement*> TracksvsEtaP;
0061   std::vector<MonitorElement*> HitsvsP;
0062   std::vector<MonitorElement*> HitsvsEta;
0063   std::vector<MonitorElement*> LayersvsP;
0064   std::vector<MonitorElement*> LayersvsEta;
0065
0066   std::vector<MonitorElement*> Num;
0067
0068   int totalNEvt;
0069
0070   const TrackerGeometry* theGeometry;
0071
0072   int numfast;
0073   int numfull;
0074   int numfastHP;
0075   int numfullHP;
0076
0077   reco::TrackBase::TrackQuality _trackQuality;
0078 };
0079
0080 testGeneralTracks::testGeneralTracks(const edm::ParameterSet& p)
0081     : geomToken(esConsumes()),
0082       pdtToken(esConsumes()),
0083       tokSimTrack_(consumes<std::vector<SimTrack>>(edm::InputTag("fastSimProducer"))),
0084       tokSimVertex_(consumes<std::vector<SimVertex>>(edm::InputTag("fastSimProducer"))),
0085       mySimEvent(2, static_cast<FSimEvent*>(0)),
0086       h0(2, static_cast<MonitorElement*>(0)),
0087       TracksvsEtaP(2, static_cast<MonitorElement*>(0)),
0088       HitsvsP(2, static_cast<MonitorElement*>(0)),
0089       HitsvsEta(2, static_cast<MonitorElement*>(0)),
0090       LayersvsP(2, static_cast<MonitorElement*>(0)),
0091       LayersvsEta(2, static_cast<MonitorElement*>(0)),
0092       Num(2, static_cast<MonitorElement*>(0)),
0093       totalNEvt(0) {
0094   // Let's just initialize the SimEvent's
0095   particleFilter_ = p.getParameter<edm::ParameterSet>("TestParticleFilter");
0096
0097   allTracks_.emplace_back(p.getParameter<edm::InputTag>("Full"));
0098   allTracks_.emplace_back(p.getParameter<edm::InputTag>("Fast"));
0099   for (const auto& tag : allTracks_)
0100     allTrackToken_.emplace_back(consumes<reco::TrackCollection>(tag));
0101
0102   // For the full sim
0103   mySimEvent[0] = new FSimEvent(particleFilter_);
0104   // For the fast sim
0105   mySimEvent[1] = new FSimEvent(particleFilter_);
0106
0107   numfast = numfull = 0;
0108   numfastHP = numfullHP = 0;
0109   _trackQuality = reco::TrackBase::qualityByName("highPurity");
0110 }
0111
0112 void testGeneralTracks::bookHistograms(DQMStore::IBooker& ibooker,
0113                                        edm::Run const& iRun,
0114                                        edm::EventSetup const& iSetup) {
0115   ibooker.setCurrentFolder("testGeneralTracks");
0116
0117   // ... and the histograms
0118   h0[0] = ibooker.book1D("generatedEta", "Generated Eta", 300, -3., 3.);
0119   h0[1] = ibooker.book1D("generatedMom", "Generated momentum", 100, 0., 10.);
0120
0121   genTracksvsEtaP = ibooker.book2D("genEtaP", "Generated eta vs p", 28, -2.8, 2.8, 100, 0, 10.);
0122   TracksvsEtaP[0] = ibooker.book2D("eff0Full", "Efficiency 0th Full", 28, -2.8, 2.8, 100, 0, 10.);
0123   TracksvsEtaP[1] = ibooker.book2D("eff0Fast", "Efficiency 0th Fast", 28, -2.8, 2.8, 100, 0, 10.);
0124   HitsvsP[0] = ibooker.book2D("Hits0PFull", "Hits vs P 0th Full", 100, 0., 10., 30, 0, 30.);
0125   HitsvsP[1] = ibooker.book2D("Hits0PFast", "Hits vs P 0th Fast", 100, 0., 10., 30, 0, 30.);
0126   HitsvsEta[0] = ibooker.book2D("Hits0EtaFull", "Hits vs Eta 0th Full", 28, -2.8, 2.8, 30, 0, 30.);
0127   HitsvsEta[1] = ibooker.book2D("Hits0EtaFast", "Hits vs Eta 0th Fast", 28, -2.8, 2.8, 30, 0, 30.);
0128   LayersvsP[0] = ibooker.book2D("Layers0PFull", "Layers vs P 0th Full", 100, 0., 10., 30, 0, 30.);
0129   LayersvsP[1] = ibooker.book2D("Layers0PFast", "Layers vs P 0th Fast", 100, 0., 10., 30, 0, 30.);
0130   LayersvsEta[0] = ibooker.book2D("Layers0EtaFull", "Layers vs Eta 0th Full", 28, -2.8, 2.8, 30, 0, 30.);
0131   LayersvsEta[1] = ibooker.book2D("Layers0EtaFast", "Layers vs Eta 0th Fast", 28, -2.8, 2.8, 30, 0, 30.);
0132 }
0133
0134 testGeneralTracks::~testGeneralTracks() {
0135   edm::LogPrint("testGeneralTracks") << "\t\t Number of Tracks ";
0136   edm::LogPrint("testGeneralTracks") << "\tFULL\t" << numfull << "\t HP= ";
0137   edm::LogPrint("testGeneralTracks") << "\tFAST\t" << numfast << "\t HP= ";
0138 }
0139
0140 void testGeneralTracks::dqmBeginRun(edm::Run const&, edm::EventSetup const& es) {
0141   // init Particle data table (from Pythia)
0142   const HepPDT::ParticleDataTable* pdt = &es.getData(pdtToken);
0143
0144   mySimEvent[0]->initializePdt(pdt);
0145   mySimEvent[1]->initializePdt(pdt);
0146
0147   theGeometry = &es.getData(geomToken);
0148 }
0149
0150 void testGeneralTracks::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0151   ++totalNEvt;
0152
0153   //  edm::LogPrint("testGeneralTracks") << " >>>>>>>>> Analizying Event " << totalNEvt << "<<<<<<< ";
0154
0155   if (totalNEvt / 1000 * 1000 == totalNEvt)
0156     edm::LogPrint("testGeneralTracks") << "Number of event analysed " << totalNEvt;
0157
0158   std::unique_ptr<edm::SimTrackContainer> nuclSimTracks(new edm::SimTrackContainer);
0159
0160   const edm::Handle<std::vector<SimTrack>>& fastSimTracks = iEvent.getHandle(tokSimTrack_);
0161   const edm::Handle<std::vector<SimVertex>>& fastSimVertices = iEvent.getHandle(tokSimVertex_);
0162   mySimEvent[1]->fill(*fastSimTracks, *fastSimVertices);
0163
0164   if (!mySimEvent[1]->nVertices())
0165     return;
0166   if (!mySimEvent[1]->nTracks())
0167     return;
0168   const FSimTrack& thePion = mySimEvent[1]->track(0);
0169
0170   double etaGen = thePion.momentum().Eta();
0171   double pGen = std::sqrt(thePion.momentum().Vect().Perp2());
0172   if (pGen < 0.2)
0173     return;
0174
0175   h0[0]->Fill(pGen);
0176   h0[1]->Fill(etaGen);
0177   genTracksvsEtaP->Fill(etaGen, pGen, 1.);
0178   //  edm::LogPrint("testGeneralTracks") << " PArticle list: Pt = "  <<  pGen << " , eta = " << etaGen << std::endl;
0179
0180   std::vector<bool> firstSeed(2, static_cast<bool>(false));
0181   std::vector<bool> secondSeed(2, static_cast<bool>(false));
0182
0183   for (unsigned ievt = 0; ievt < allTracks_.size(); ++ievt) {
0184     const edm::Handle<reco::TrackCollection>& tkRef0 = iEvent.getHandle(allTrackToken_[ievt]);
0185     std::vector<const reco::TrackCollection*> tkColl;
0186     tkColl.push_back(tkRef0.product());
0187
0188     reco::TrackCollection::const_iterator itk0 = tkColl[0]->begin();
0189     reco::TrackCollection::const_iterator itk0_e = tkColl[0]->end();
0190     for (; itk0 != itk0_e; ++itk0) {
0191       //edm::LogPrint("testGeneralTracks") << "quality " << itk0->quality(_trackQuality) << std::endl;
0192       if (!(itk0->quality(_trackQuality))) {
0193         //edm::LogPrint("testGeneralTracks") << "evt " << totalNEvt << "\tTRACK REMOVED" << std::endl;
0194         continue;
0195       }
0196       if (ievt == 0)
0197         numfullHP++;
0198       if (ievt == 1)
0199         numfastHP++;
0200       TracksvsEtaP[ievt]->Fill(etaGen, pGen, 1.);
0201       HitsvsEta[ievt]->Fill(etaGen, itk0->found(), 1.);
0202       HitsvsP[ievt]->Fill(pGen, itk0->found(), 1.);
0203       LayersvsEta[ievt]->Fill(etaGen, itk0->hitPattern().trackerLayersWithMeasurement(), 1.);
0204       LayersvsP[ievt]->Fill(pGen, itk0->hitPattern().trackerLayersWithMeasurement(), 1.);
0205     }
0206
0207     //    edm::LogPrint("testGeneralTracks") << "\t\t Number of Tracks " << std::endl;
0208     if (ievt == 0) {
0209       numfull += tkColl[0]->size();
0210       //      edm::LogPrint("testGeneralTracks") << "\tFULL\t" << tkColl[0]->size() << "\t" << numfull << "\t" << numfullHP << std::endl;
0211     } else if (ievt == 1) {
0212       numfast += tkColl[0]->size();
0213       // edm::LogPrint("testGeneralTracks") << "\tFAST\t" << tkColl[0]->size() << "\t" << numfast << "\t" << numfastHP << std::endl;
0214     }
0215   }
0216 }
0217
0218 //define this as a plug-in
0219
0220 DEFINE_FWK_MODULE(testGeneralTracks);