File indexing completed on 2024-09-10 02:59:11
0001
0002
0003
0004
0005
0006
0007 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0008 #include "SimMuon/DTDigitizer/test/DTDigiAnalyzer.h"
0009 #include <DataFormats/DTDigi/interface/DTDigiCollection.h>
0010 #include <FWCore/Framework/interface/Event.h>
0011
0012
0013
0014
0015
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020
0021 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0022 #include "Geometry/DTGeometry/interface/DTLayer.h"
0023
0024 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0025 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0026
0027 #include <iostream>
0028 #include <string>
0029
0030 #include "TFile.h"
0031
0032 #include "SimMuon/DTDigitizer/interface/Histograms.h"
0033
0034 using namespace edm;
0035 using namespace std;
0036
0037 DTDigiAnalyzer::DTDigiAnalyzer(const ParameterSet &pset)
0038 : hDigis_global("Global"), hDigis_W0("Wheel0"), hDigis_W1("Wheel1"), hDigis_W2("Wheel2"), hAllHits("AllHits") {
0039
0040
0041
0042 label = pset.getUntrackedParameter<string>("label");
0043 tokGeo_ = esConsumes<DTGeometry, MuonGeometryRecord>();
0044 file = new TFile("DTDigiPlots.root", "RECREATE");
0045 file->cd();
0046 DigiTimeBox = new TH1F("DigiTimeBox", "Digi Time Box", 2048, 0, 1600);
0047 if (file->IsOpen())
0048 cout << "file open!" << endl;
0049 else
0050 cout << "*** Error in opening file ***" << endl;
0051
0052 psim_token = consumes<PSimHitContainer>(edm::InputTag("g4SimHits", "MuonDTHits"));
0053 DTd_token = consumes<DTDigiCollection>(edm::InputTag(label));
0054 }
0055
0056 DTDigiAnalyzer::~DTDigiAnalyzer() {}
0057
0058 void DTDigiAnalyzer::endJob() {
0059
0060
0061 file->cd();
0062 DigiTimeBox->Write();
0063 hDigis_global.Write();
0064 hDigis_W0.Write();
0065 hDigis_W1.Write();
0066 hDigis_W2.Write();
0067 hAllHits.Write();
0068 file->Close();
0069
0070
0071 }
0072
0073 void DTDigiAnalyzer::analyze(const Event &event, const EventSetup &eventSetup) {
0074 cout << "--- Run: " << event.id().run() << " Event: " << event.id().event() << endl;
0075
0076 Handle<DTDigiCollection> dtDigis;
0077 event.getByToken(DTd_token, dtDigis);
0078
0079 Handle<PSimHitContainer> simHits;
0080 event.getByToken(psim_token, simHits);
0081
0082 auto muonGeom = eventSetup.getHandle(tokGeo_);
0083
0084 DTWireIdMap wireMap;
0085
0086 for (vector<PSimHit>::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++) {
0087
0088 DTWireId wireId(hit->detUnitId());
0089
0090 wireMap[wireId].push_back(&(*hit));
0091
0092 LocalPoint entryP = hit->entryPoint();
0093 LocalPoint exitP = hit->exitPoint();
0094 int partType = hit->particleType();
0095
0096 float path = (exitP - entryP).mag();
0097 float path_x = fabs((exitP - entryP).x());
0098
0099 hAllHits.Fill(entryP.x(),
0100 exitP.x(),
0101 entryP.y(),
0102 exitP.y(),
0103 entryP.z(),
0104 exitP.z(),
0105 path,
0106 path_x,
0107 partType,
0108 hit->processType(),
0109 hit->pabs());
0110
0111 if (hit->timeOfFlight() > 1e4) {
0112 cout << "PID: " << hit->particleType() << " TOF: " << hit->timeOfFlight() << " Proc Type: " << hit->processType()
0113 << " p: " << hit->pabs() << endl;
0114 hAllHits.FillTOF(hit->timeOfFlight());
0115 }
0116 }
0117
0118 DTDigiCollection::DigiRangeIterator detUnitIt;
0119 for (detUnitIt = dtDigis->begin(); detUnitIt != dtDigis->end(); ++detUnitIt) {
0120 const DTLayerId &id = (*detUnitIt).first;
0121 const DTDigiCollection::Range &range = (*detUnitIt).second;
0122
0123
0124 cout << "--------------" << endl;
0125 cout << "id: " << id;
0126
0127
0128 for (DTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0129 cout << " Wire: " << (*digiIt).wire() << endl << " digi time (ns): " << (*digiIt).time() << endl;
0130 DigiTimeBox->Fill((*digiIt).time());
0131
0132 DTWireId wireId(id, (*digiIt).wire());
0133 int mu = 0;
0134 float theta = 0;
0135
0136 for (vector<const PSimHit *>::iterator hit = wireMap[wireId].begin(); hit != wireMap[wireId].end(); hit++) {
0137 cout << "momentum x: " << (*hit)->momentumAtEntry().x() << endl
0138 << "momentum z: " << (*hit)->momentumAtEntry().z() << endl;
0139 if (abs((*hit)->particleType()) == 13) {
0140 theta = atan((*hit)->momentumAtEntry().x() / (-(*hit)->momentumAtEntry().z())) * 180 / M_PI;
0141 cout << "atan: " << theta << endl;
0142 mu++;
0143 } else {
0144
0145
0146
0147
0148 }
0149 }
0150 if (mu && theta) {
0151 hDigis_global.Fill((*digiIt).time(), theta, id.superlayer());
0152
0153 WheelHistos(id.wheel())->Fill((*digiIt).time(), theta, id.superlayer());
0154 }
0155
0156 }
0157 }
0158 cout << "--------------" << endl;
0159 }
0160
0161 hDigis *DTDigiAnalyzer::WheelHistos(int wheel) {
0162 switch (abs(wheel)) {
0163 case 0:
0164 return &hDigis_W0;
0165
0166 case 1:
0167 return &hDigis_W1;
0168
0169 case 2:
0170 return &hDigis_W2;
0171
0172 default:
0173 return nullptr;
0174 }
0175 }
0176
0177 #include "FWCore/Framework/interface/MakerMacros.h"
0178 #include "FWCore/PluginManager/interface/ModuleDef.h"
0179
0180 DEFINE_FWK_MODULE(DTDigiAnalyzer);