DTDigitizer/test/DTDigiAnalyzer.cc

0001 /** \class DTDigiAnalyzer
0002  *  Analyse the the muon-drift-tubes digitizer.
0003  *
0004  *  \authors: R. Bellan
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 // #include "SimMuon/DTDigitizer/test/analysis/DTMCStatistics.h"
0013 // #include "SimMuon/DTDigitizer/test/analysis/DTMuonDigiStatistics.h"
0014 // #include "SimMuon/DTDigitizer/test/analysis/DTHitsAnalysis.h"
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   //   MCStatistics = new DTMCStatistics();
0040   //   MuonDigiStatistics = new DTMuonDigiStatistics();
0041   //   HitsAnalysis = new DTHitsAnalysis();
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   // cout<<"Number of analyzed event: "<<nevts<<endl;
0060   // HitsAnalysis->Report();
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   //    delete file;
0070   // delete DigiTimeBox;
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     // Create the id of the wire, the simHits in the DT known also the wireId
0088     DTWireId wireId(hit->detUnitId());
0089     // Fill the map
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     // DTLayerId print-out
0124     cout << "--------------" << endl;
0125     cout << "id: " << id;
0126
0127     // Loop over the digis of this DetUnit
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           //      cout<<"PID: "<<(*hit)->particleType()
0145           //          <<" TOF: "<<(*hit)->timeOfFlight()
0146           //          <<" Proc Type: "<<(*hit)->processType()
0147           //          <<" p: " << (*hit)->pabs() <<endl;
0148         }
0149       }
0150       if (mu && theta) {
0151         hDigis_global.Fill((*digiIt).time(), theta, id.superlayer());
0152         // filling digi histos for wheel and for RZ and RPhi
0153         WheelHistos(id.wheel())->Fill((*digiIt).time(), theta, id.superlayer());
0154       }
0155
0156     }  // for digis in layer
0157   }  // for layers
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);