DTDigitizer/src/DTDigiReader.cc

0001 #ifndef SimMuon_DTDigiReader_h
0002 #define SimMuon_DTDigiReader_h
0003
0004 /** \class DTDigiReader
0005  *  Analyse the the muon-drift-tubes digitizer.
0006  *
0007  *  \authors: R. Bellan
0008  */
0009
0010 #include "FWCore/Framework/interface/ConsumesCollector.h"
0011
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include <FWCore/Framework/interface/Event.h>
0014 #include <FWCore/Framework/interface/one/EDAnalyzer.h>
0015
0016 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0017 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0018
0019 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0020 #include <DataFormats/DTDigi/interface/DTDigiCollection.h>
0021
0022 #include <iostream>
0023
0024 #include "TFile.h"
0025 #include "TH1F.h"  //FIXME
0026
0027 using namespace edm;
0028 using namespace std;
0029
0030 class DTDigiReader : public edm::one::EDAnalyzer<> {
0031 public:
0032   explicit DTDigiReader(const ParameterSet &pset) {
0033     file = new TFile("DTDigiPlots.root", "RECREATE");
0034     file->cd();
0035     DigiTimeBox = new TH1F("DigiTimeBox", "Digi Time Box", 2048, 0, 1600);
0036     DigiTimeBoxW0 = new TH1F("DigiTimeBoxW0", "Digi Time Box W0", 2000, 0, 1600);
0037     DigiTimeBoxW1 = new TH1F("DigiTimeBoxW1", "Digi Time Box W1", 2000, 0, 1600);
0038     DigiTimeBoxW2 = new TH1F("DigiTimeBoxW2", "Digi Time Box W2", 2000, 0, 1600);
0039     if (file->IsOpen())
0040       cout << "file open!" << endl;
0041     else
0042       cout << "*** Error in opening file ***" << endl;
0043     label = pset.getUntrackedParameter<string>("label");
0044     psim_token = consumes<PSimHitContainer>(edm::InputTag("g4SimHits", "MuonDTHits"));
0045     DTd_token = consumes<DTDigiCollection>(edm::InputTag(label));
0046   }
0047
0048   ~DTDigiReader() override {
0049     file->cd();
0050     DigiTimeBox->Write();
0051     DigiTimeBoxW0->Write();
0052     DigiTimeBoxW1->Write();
0053     DigiTimeBoxW2->Write();
0054     file->Close();
0055     //    delete file;
0056     // delete DigiTimeBox;
0057   }
0058
0059   void analyze(const Event &event, const EventSetup &eventSetup) override {
0060     cout << "--- Run: " << event.id().run() << " Event: " << event.id().event() << endl;
0061
0062     Handle<DTDigiCollection> dtDigis;
0063     event.getByToken(DTd_token, dtDigis);
0064     // event.getByLabel("MuonDTDigis", dtDigis);
0065     Handle<PSimHitContainer> simHits;
0066     event.getByToken(psim_token, simHits);
0067
0068     DTDigiCollection::DigiRangeIterator detUnitIt;
0069     for (detUnitIt = dtDigis->begin(); detUnitIt != dtDigis->end(); ++detUnitIt) {
0070       const DTLayerId &id = (*detUnitIt).first;
0071       const DTDigiCollection::Range &range = (*detUnitIt).second;
0072
0073       // DTLayerId print-out
0074       cout << "--------------" << endl;
0075       cout << "id: " << id;
0076
0077       // Loop over the digis of this DetUnit
0078       for (DTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0079         //  if((*digiIt).time()<703 &&(*digiIt).time()>699) {
0080         cout << " Wire: " << (*digiIt).wire() << endl << " digi time (ns): " << (*digiIt).time() << endl;
0081
0082         for (vector<PSimHit>::const_iterator simHit = simHits->begin(); simHit != simHits->end(); simHit++) {
0083           DTWireId wireId((*simHit).detUnitId());
0084           if (wireId.layerId() == id && abs((*simHit).particleType()) == 13) {
0085             cout << "entry: " << (*simHit).entryPoint() << endl
0086                  << "exit: " << (*simHit).exitPoint() << endl
0087                  << "TOF: " << (*simHit).timeOfFlight() << endl;
0088           }
0089         }
0090
0091         //  }
0092
0093         if (id.layer() == 3)
0094           DigiTimeBoxW0->Fill((*digiIt).time());
0095         else if (abs(id.superlayer()) == 1)
0096           DigiTimeBoxW1->Fill((*digiIt).time());
0097         else if (abs(id.superlayer()) == 2)
0098           DigiTimeBoxW2->Fill((*digiIt).time());
0099         else
0100           cout << "Error" << endl;
0101         DigiTimeBox->Fill((*digiIt).time());
0102
0103       }  // for digis in layer
0104     }    // for layers
0105     cout << "--------------" << endl;
0106   }
0107
0108 private:
0109   string label;
0110   TH1F *DigiTimeBox;
0111   TH1F *DigiTimeBoxW0;
0112   TH1F *DigiTimeBoxW1;
0113   TH1F *DigiTimeBoxW2;
0114   TFile *file;
0115
0116   edm::EDGetTokenT<PSimHitContainer> psim_token;
0117   edm::EDGetTokenT<DTDigiCollection> DTd_token;
0118 };
0119
0120 #endif