Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-10 23:39:21

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   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken;
0046   edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord> pdtToken;
0047 
0048   // See RecoParticleFlow/PFProducer/interface/PFProducer.h
0049   edm::ParameterSet particleFilter_;
0050   std::vector<edm::InputTag> allTracks;
0051   bool saveNU;
0052   std::vector<FSimEvent*> mySimEvent;
0053   std::string simModuleLabel_;
0054   // Histograms
0055   std::vector<MonitorElement*> h0;
0056   MonitorElement* genTracksvsEtaP;
0057   std::vector<MonitorElement*> TracksvsEtaP;
0058   std::vector<MonitorElement*> HitsvsP;
0059   std::vector<MonitorElement*> HitsvsEta;
0060   std::vector<MonitorElement*> LayersvsP;
0061   std::vector<MonitorElement*> LayersvsEta;
0062 
0063   std::vector<MonitorElement*> Num;
0064 
0065   int totalNEvt;
0066 
0067   const TrackerGeometry* theGeometry;
0068 
0069   int numfast;
0070   int numfull;
0071   int numfastHP;
0072   int numfullHP;
0073 
0074   reco::TrackBase::TrackQuality _trackQuality;
0075 };
0076 
0077 testGeneralTracks::testGeneralTracks(const edm::ParameterSet& p)
0078     : geomToken(esConsumes()),
0079       pdtToken(esConsumes()),
0080       mySimEvent(2, static_cast<FSimEvent*>(0)),
0081       h0(2, static_cast<MonitorElement*>(0)),
0082       TracksvsEtaP(2, static_cast<MonitorElement*>(0)),
0083       HitsvsP(2, static_cast<MonitorElement*>(0)),
0084       HitsvsEta(2, static_cast<MonitorElement*>(0)),
0085       LayersvsP(2, static_cast<MonitorElement*>(0)),
0086       LayersvsEta(2, static_cast<MonitorElement*>(0)),
0087 
0088       Num(2, static_cast<MonitorElement*>(0)),
0089 
0090       totalNEvt(0) {
0091   // Let's just initialize the SimEvent's
0092   particleFilter_ = p.getParameter<edm::ParameterSet>("TestParticleFilter");
0093 
0094   allTracks.push_back(p.getParameter<edm::InputTag>("Full"));
0095   allTracks.push_back(p.getParameter<edm::InputTag>("Fast"));
0096 
0097   // For the full sim
0098   mySimEvent[0] = new FSimEvent(particleFilter_);
0099   // For the fast sim
0100   mySimEvent[1] = new FSimEvent(particleFilter_);
0101 
0102   numfast = numfull = 0;
0103   numfastHP = numfullHP = 0;
0104   _trackQuality = reco::TrackBase::qualityByName("highPurity");
0105 }
0106 
0107 void testGeneralTracks::bookHistograms(DQMStore::IBooker& ibooker,
0108                                        edm::Run const& iRun,
0109                                        edm::EventSetup const& iSetup) {
0110   ibooker.setCurrentFolder("testGeneralTracks");
0111 
0112   // ... and the histograms
0113   h0[0] = ibooker.book1D("generatedEta", "Generated Eta", 300, -3., 3.);
0114   h0[1] = ibooker.book1D("generatedMom", "Generated momentum", 100, 0., 10.);
0115 
0116   genTracksvsEtaP = ibooker.book2D("genEtaP", "Generated eta vs p", 28, -2.8, 2.8, 100, 0, 10.);
0117   TracksvsEtaP[0] = ibooker.book2D("eff0Full", "Efficiency 0th Full", 28, -2.8, 2.8, 100, 0, 10.);
0118   TracksvsEtaP[1] = ibooker.book2D("eff0Fast", "Efficiency 0th Fast", 28, -2.8, 2.8, 100, 0, 10.);
0119   HitsvsP[0] = ibooker.book2D("Hits0PFull", "Hits vs P 0th Full", 100, 0., 10., 30, 0, 30.);
0120   HitsvsP[1] = ibooker.book2D("Hits0PFast", "Hits vs P 0th Fast", 100, 0., 10., 30, 0, 30.);
0121   HitsvsEta[0] = ibooker.book2D("Hits0EtaFull", "Hits vs Eta 0th Full", 28, -2.8, 2.8, 30, 0, 30.);
0122   HitsvsEta[1] = ibooker.book2D("Hits0EtaFast", "Hits vs Eta 0th Fast", 28, -2.8, 2.8, 30, 0, 30.);
0123   LayersvsP[0] = ibooker.book2D("Layers0PFull", "Layers vs P 0th Full", 100, 0., 10., 30, 0, 30.);
0124   LayersvsP[1] = ibooker.book2D("Layers0PFast", "Layers vs P 0th Fast", 100, 0., 10., 30, 0, 30.);
0125   LayersvsEta[0] = ibooker.book2D("Layers0EtaFull", "Layers vs Eta 0th Full", 28, -2.8, 2.8, 30, 0, 30.);
0126   LayersvsEta[1] = ibooker.book2D("Layers0EtaFast", "Layers vs Eta 0th Fast", 28, -2.8, 2.8, 30, 0, 30.);
0127 }
0128 
0129 testGeneralTracks::~testGeneralTracks() {
0130   edm::LogPrint("testGeneralTracks") << "\t\t Number of Tracks ";
0131   edm::LogPrint("testGeneralTracks") << "\tFULL\t" << numfull << "\t HP= ";
0132   edm::LogPrint("testGeneralTracks") << "\tFAST\t" << numfast << "\t HP= ";
0133 }
0134 
0135 void testGeneralTracks::dqmBeginRun(edm::Run const&, edm::EventSetup const& es) {
0136   // init Particle data table (from Pythia)
0137   const HepPDT::ParticleDataTable* pdt = &es.getData(pdtToken);
0138 
0139   mySimEvent[0]->initializePdt(pdt);
0140   mySimEvent[1]->initializePdt(pdt);
0141 
0142   theGeometry = &es.getData(geomToken);
0143 }
0144 
0145 void testGeneralTracks::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0146   ++totalNEvt;
0147 
0148   //  edm::LogPrint("testGeneralTracks") << " >>>>>>>>> Analizying Event " << totalNEvt << "<<<<<<< ";
0149 
0150   if (totalNEvt / 1000 * 1000 == totalNEvt)
0151     edm::LogPrint("testGeneralTracks") << "Number of event analysed " << totalNEvt;
0152 
0153   std::unique_ptr<edm::SimTrackContainer> nuclSimTracks(new edm::SimTrackContainer);
0154 
0155   edm::Handle<std::vector<SimTrack> > fastSimTracks;
0156   iEvent.getByLabel("fastSimProducer", fastSimTracks);
0157   edm::Handle<std::vector<SimVertex> > fastSimVertices;
0158   iEvent.getByLabel("fastSimProducer", fastSimVertices);
0159   mySimEvent[1]->fill(*fastSimTracks, *fastSimVertices);
0160 
0161   if (!mySimEvent[1]->nVertices())
0162     return;
0163   if (!mySimEvent[1]->nTracks())
0164     return;
0165   const FSimTrack& thePion = mySimEvent[1]->track(0);
0166 
0167   double etaGen = thePion.momentum().Eta();
0168   double pGen = std::sqrt(thePion.momentum().Vect().Perp2());
0169   if (pGen < 0.2)
0170     return;
0171 
0172   h0[0]->Fill(pGen);
0173   h0[1]->Fill(etaGen);
0174   genTracksvsEtaP->Fill(etaGen, pGen, 1.);
0175   //  edm::LogPrint("testGeneralTracks") << " PArticle list: Pt = "  <<  pGen << " , eta = " << etaGen << std::endl;
0176 
0177   std::vector<bool> firstSeed(2, static_cast<bool>(false));
0178   std::vector<bool> secondSeed(2, static_cast<bool>(false));
0179 
0180   for (unsigned ievt = 0; ievt < 2; ++ievt) {
0181     edm::Handle<reco::TrackCollection> tkRef0;
0182     iEvent.getByLabel(allTracks[ievt], tkRef0);
0183     std::vector<const reco::TrackCollection*> tkColl;
0184     tkColl.push_back(tkRef0.product());
0185 
0186     reco::TrackCollection::const_iterator itk0 = tkColl[0]->begin();
0187     reco::TrackCollection::const_iterator itk0_e = tkColl[0]->end();
0188     for (; itk0 != itk0_e; ++itk0) {
0189       //edm::LogPrint("testGeneralTracks") << "quality " << itk0->quality(_trackQuality) << std::endl;
0190       if (!(itk0->quality(_trackQuality))) {
0191         //edm::LogPrint("testGeneralTracks") << "evt " << totalNEvt << "\tTRACK REMOVED" << std::endl;
0192         continue;
0193       }
0194       if (ievt == 0)
0195         numfullHP++;
0196       if (ievt == 1)
0197         numfastHP++;
0198       TracksvsEtaP[ievt]->Fill(etaGen, pGen, 1.);
0199       HitsvsEta[ievt]->Fill(etaGen, itk0->found(), 1.);
0200       HitsvsP[ievt]->Fill(pGen, itk0->found(), 1.);
0201       LayersvsEta[ievt]->Fill(etaGen, itk0->hitPattern().trackerLayersWithMeasurement(), 1.);
0202       LayersvsP[ievt]->Fill(pGen, itk0->hitPattern().trackerLayersWithMeasurement(), 1.);
0203     }
0204 
0205     //    edm::LogPrint("testGeneralTracks") << "\t\t Number of Tracks " << std::endl;
0206     if (ievt == 0) {
0207       numfull += tkColl[0]->size();
0208       //      edm::LogPrint("testGeneralTracks") << "\tFULL\t" << tkColl[0]->size() << "\t" << numfull << "\t" << numfullHP << std::endl;
0209     } else if (ievt == 1) {
0210       numfast += tkColl[0]->size();
0211       // edm::LogPrint("testGeneralTracks") << "\tFAST\t" << tkColl[0]->size() << "\t" << numfast << "\t" << numfastHP << std::endl;
0212     }
0213   }
0214 }
0215 
0216 //define this as a plug-in
0217 
0218 DEFINE_FWK_MODULE(testGeneralTracks);