CMS 3D CMS Logo

MuonDTDigis.cc
Go to the documentation of this file.
1 
7 #include "MuonDTDigis.h"
8 
9 #include <iostream>
10 #include <string>
11 
12 #include "TFile.h"
13 
14 #include "SimMuon/DTDigitizer/test/Histograms.h"
15 
16 using namespace edm;
17 using namespace std;
18 
20  // ----------------------
21  // Get the debug parameter for verbose output
22  verbose_ = pset.getUntrackedParameter<bool>("verbose", false);
23 
24  // the name of the Digi collection
25  SimHitToken_ = consumes<edm::PSimHitContainer>(pset.getParameter<edm::InputTag>("SimHitLabel"));
26 
27  // the name of the Digi collection
28  DigiToken_ = consumes<DTDigiCollection>(pset.getParameter<edm::InputTag>(("DigiLabel")));
29 
30  hDigis_global = std::make_unique<hDigis>("Global");
31  hDigis_W0 = std::make_unique<hDigis>("Wheel0");
32  hDigis_W1 = std::make_unique<hDigis>("Wheel1");
33  hDigis_W2 = std::make_unique<hDigis>("Wheel2");
34  hAllHits = std::make_unique<hHits>("AllHits");
35 }
36 
38  if (verbose_)
39  cout << "[MuonDTDigis] Destructor called" << endl;
40 }
41 
43  edm::Run const &iRun,
44  edm::EventSetup const & /* iSetup */) {
45  meDigiTimeBox_ = nullptr;
46  meDigiTimeBox_wheel2m_ = nullptr;
47  meDigiTimeBox_wheel1m_ = nullptr;
48  meDigiTimeBox_wheel0_ = nullptr;
49  meDigiTimeBox_wheel1p_ = nullptr;
50  meDigiTimeBox_wheel2p_ = nullptr;
51  meDigiEfficiency_ = nullptr;
52  meDigiEfficiencyMu_ = nullptr;
53  meDoubleDigi_ = nullptr;
54  meSimvsDigi_ = nullptr;
55  meWire_DoubleDigi_ = nullptr;
56 
57  meMB1_sim_occup_ = nullptr;
58  meMB1_digi_occup_ = nullptr;
59  meMB2_sim_occup_ = nullptr;
60  meMB2_digi_occup_ = nullptr;
61  meMB3_sim_occup_ = nullptr;
62  meMB3_digi_occup_ = nullptr;
63  meMB4_sim_occup_ = nullptr;
64  meMB4_digi_occup_ = nullptr;
65 
66  meDigiHisto_ = nullptr;
67 
68  // ----------------------
69  // We go
70  // ----------------------
71 
72  Char_t histo_n[100];
73  Char_t histo_t[100];
74 
75  iBooker.setCurrentFolder("MuonDTDigisV/DTDigiValidationTask");
76 
77  sprintf(histo_n, "DigiTimeBox");
78  sprintf(histo_t, "Digi Time Box");
79  meDigiTimeBox_ = iBooker.book1D(histo_n, histo_t, 1536, 0, 1200);
80 
81  sprintf(histo_n, "DigiTimeBox_wheel2m");
82  sprintf(histo_t, "Digi Time Box wheel -2");
83  meDigiTimeBox_wheel2m_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
84 
85  sprintf(histo_n, "DigiTimeBox_wheel1m");
86  sprintf(histo_t, "Digi Time Box wheel -1");
87  meDigiTimeBox_wheel1m_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
88 
89  sprintf(histo_n, "DigiTimeBox_wheel0");
90  sprintf(histo_t, "Digi Time Box wheel 0");
91  meDigiTimeBox_wheel0_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
92 
93  sprintf(histo_n, "DigiTimeBox_wheel1p");
94  sprintf(histo_t, "Digi Time Box wheel 1");
95  meDigiTimeBox_wheel1p_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
96 
97  sprintf(histo_n, "DigiTimeBox_wheel2p");
98  sprintf(histo_t, "Digi Time Box wheel 2");
99  meDigiTimeBox_wheel2p_ = iBooker.book1D(histo_n, histo_t, 384, 0, 1200);
100 
101  sprintf(histo_n, "DigiEfficiencyMu");
102  sprintf(histo_t, "Ratio (#Digis Mu)/(#SimHits Mu)");
103  meDigiEfficiencyMu_ = iBooker.book1D(histo_n, histo_t, 100, 0., 5.);
104 
105  sprintf(histo_n, "DigiEfficiency");
106  sprintf(histo_t, "Ratio (#Digis)/(#SimHits)");
107  meDigiEfficiency_ = iBooker.book1D(histo_n, histo_t, 100, 0., 5.);
108 
109  sprintf(histo_n, "Number_Digi_per_layer");
110  sprintf(histo_t, "Number_Digi_per_layer");
111  meDoubleDigi_ = iBooker.book1D(histo_n, histo_t, 10, 0., 10.);
112 
113  sprintf(histo_n, "Number_simhit_vs_digi");
114  sprintf(histo_t, "Number_simhit_vs_digi");
115  meSimvsDigi_ = iBooker.book2D(histo_n, histo_t, 100, 0., 140., 100, 0., 140.);
116 
117  sprintf(histo_n, "Wire_Number_with_double_Digi");
118  sprintf(histo_t, "Wire_Number_with_double_Digi");
119  meWire_DoubleDigi_ = iBooker.book1D(histo_n, histo_t, 100, 0., 100.);
120 
121  sprintf(histo_n, "Simhit_occupancy_MB1");
122  sprintf(histo_t, "Simhit_occupancy_MB1");
123  meMB1_sim_occup_ = iBooker.book1D(histo_n, histo_t, 55, 0., 55.);
124 
125  sprintf(histo_n, "Digi_occupancy_MB1");
126  sprintf(histo_t, "Digi_occupancy_MB1");
127  meMB1_digi_occup_ = iBooker.book1D(histo_n, histo_t, 55, 0., 55.);
128 
129  sprintf(histo_n, "Simhit_occupancy_MB2");
130  sprintf(histo_t, "Simhit_occupancy_MB2");
131  meMB2_sim_occup_ = iBooker.book1D(histo_n, histo_t, 63, 0., 63.);
132 
133  sprintf(histo_n, "Digi_occupancy_MB2");
134  sprintf(histo_t, "Digi_occupancy_MB2");
135  meMB2_digi_occup_ = iBooker.book1D(histo_n, histo_t, 63, 0., 63.);
136 
137  sprintf(histo_n, "Simhit_occupancy_MB3");
138  sprintf(histo_t, "Simhit_occupancy_MB3");
139  meMB3_sim_occup_ = iBooker.book1D(histo_n, histo_t, 75, 0., 75.);
140 
141  sprintf(histo_n, "Digi_occupancy_MB3");
142  sprintf(histo_t, "Digi_occupancy_MB3");
143  meMB3_digi_occup_ = iBooker.book1D(histo_n, histo_t, 75, 0., 75.);
144 
145  sprintf(histo_n, "Simhit_occupancy_MB4");
146  sprintf(histo_t, "Simhit_occupancy_MB4");
147  meMB4_sim_occup_ = iBooker.book1D(histo_n, histo_t, 99, 0., 99.);
148 
149  sprintf(histo_n, "Digi_occupancy_MB4");
150  sprintf(histo_t, "Digi_occupancy_MB4");
151  meMB4_digi_occup_ = iBooker.book1D(histo_n, histo_t, 99, 0., 99.);
152 
153  // Begona's option
154  char stringcham[40];
155  for (int slnum = 1; slnum < 62; ++slnum) {
156  sprintf(stringcham, "DigiTimeBox_slid_%d", slnum);
157  meDigiHisto_ = iBooker.book1D(stringcham, stringcham, 100, 0, 1200);
158  meDigiTimeBox_SL_.push_back(meDigiHisto_);
159  }
160 }
161 
162 void MuonDTDigis::analyze(const Event &event, const EventSetup &eventSetup) {
163  if (verbose_)
164  cout << "--- [MuonDTDigis] Analysing Event: #Run: " << event.id().run() << " #Event: " << event.id().event()
165  << endl;
166 
167  // Get the DT Geometry
168  ESHandle<DTGeometry> muonGeom;
169  eventSetup.get<MuonGeometryRecord>().get(muonGeom);
170 
171  // Get the Digi collection from the event
172  Handle<DTDigiCollection> dtDigis;
173  event.getByToken(DigiToken_, dtDigis);
174 
175  // Get simhits
177  event.getByToken(SimHitToken_, simHits);
178 
179  int num_mudigis;
180  int num_musimhits;
181  int num_digis;
182  int num_digis_layer;
183  int cham_num;
184  int wire_touched;
185  num_digis = 0;
186  num_mudigis = 0;
187  num_musimhits = 0;
188  DTWireIdMap wireMap;
189 
190  for (vector<PSimHit>::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++) {
191  // Create the id of the wire, the simHits in the DT known also the wireId
192  DTWireId wireId(hit->detUnitId());
193  // cout << " PSimHits wire id " << wireId << " part type " <<
194  // hit->particleType() << endl;
195 
196  // Fill the map
197  wireMap[wireId].push_back(&(*hit));
198 
199  LocalPoint entryP = hit->entryPoint();
200  LocalPoint exitP = hit->exitPoint();
201  int partType = hit->particleType();
202  if (abs(partType) == 13)
203  num_musimhits++;
204 
205  if (wireId.station() == 1 && abs(partType) == 13)
206  meMB1_sim_occup_->Fill(wireId.wire());
207  if (wireId.station() == 2 && abs(partType) == 13)
208  meMB2_sim_occup_->Fill(wireId.wire());
209  if (wireId.station() == 3 && abs(partType) == 13)
210  meMB3_sim_occup_->Fill(wireId.wire());
211  if (wireId.station() == 4 && abs(partType) == 13)
212  meMB4_sim_occup_->Fill(wireId.wire());
213 
214  float path = (exitP - entryP).mag();
215  float path_x = fabs((exitP - entryP).x());
216 
217  hAllHits->Fill(entryP.x(),
218  exitP.x(),
219  entryP.y(),
220  exitP.y(),
221  entryP.z(),
222  exitP.z(),
223  path,
224  path_x,
225  partType,
226  hit->processType(),
227  hit->pabs());
228  }
229 
230  // cout << "num muon simhits " << num_musimhits << endl;
231 
233  for (detUnitIt = dtDigis->begin(); detUnitIt != dtDigis->end(); ++detUnitIt) {
234  const DTLayerId &id = (*detUnitIt).first;
235  const DTDigiCollection::Range &range = (*detUnitIt).second;
236 
237  num_digis_layer = 0;
238  cham_num = 0;
239  wire_touched = 0;
240 
241  // Loop over the digis of this DetUnit
242  for (DTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
243  // cout<<" Wire: "<<(*digiIt).wire()<<endl
244  // <<" digi time (ns): "<<(*digiIt).time()<<endl;
245 
246  num_digis++;
247  num_digis_layer++;
248  if (num_digis_layer > 1) {
249  if ((*digiIt).wire() == wire_touched) {
250  meWire_DoubleDigi_->Fill((*digiIt).wire());
251  // cout << "old & new wire " << wire_touched << " " <<
252  // (*digiIt).wire() << endl;
253  }
254  }
255  wire_touched = (*digiIt).wire();
256 
257  meDigiTimeBox_->Fill((*digiIt).time());
258  if (id.wheel() == -2)
259  meDigiTimeBox_wheel2m_->Fill((*digiIt).time());
260  if (id.wheel() == -1)
261  meDigiTimeBox_wheel1m_->Fill((*digiIt).time());
262  if (id.wheel() == 0)
263  meDigiTimeBox_wheel0_->Fill((*digiIt).time());
264  if (id.wheel() == 1)
265  meDigiTimeBox_wheel1p_->Fill((*digiIt).time());
266  if (id.wheel() == 2)
267  meDigiTimeBox_wheel2p_->Fill((*digiIt).time());
268 
269  // Superlayer number and fill histo with digi timebox
270 
271  cham_num = (id.wheel() + 2) * 12 + (id.station() - 1) * 3 + id.superlayer();
272  // cout << " Histo number " << cham_num << endl;
273 
274  meDigiTimeBox_SL_[cham_num]->Fill((*digiIt).time());
275 
276  // cout << " size de digis " << (*digiIt).size() << endl;
277 
278  DTWireId wireId(id, (*digiIt).wire());
279  if (wireId.station() == 1)
280  meMB1_digi_occup_->Fill((*digiIt).wire());
281  if (wireId.station() == 2)
282  meMB2_digi_occup_->Fill((*digiIt).wire());
283  if (wireId.station() == 3)
284  meMB3_digi_occup_->Fill((*digiIt).wire());
285  if (wireId.station() == 4)
286  meMB4_digi_occup_->Fill((*digiIt).wire());
287 
288  int mu = 0;
289  float theta = 0;
290 
291  for (vector<const PSimHit *>::iterator hit = wireMap[wireId].begin(); hit != wireMap[wireId].end(); hit++)
292  if (abs((*hit)->particleType()) == 13) {
293  theta = atan((*hit)->momentumAtEntry().x() / (-(*hit)->momentumAtEntry().z())) * 180 / M_PI;
294  // cout<<"momentum x: "<<(*hit)->momentumAtEntry().x()<<endl
295  // <<"momentum z: "<<(*hit)->momentumAtEntry().z()<<endl
296  // <<"atan: "<<theta<<endl;
297  mu++;
298  }
299 
300  if (mu)
301  num_mudigis++;
302 
303  if (mu && theta) {
304  hDigis_global->Fill((*digiIt).time(), theta, id.superlayer());
305  // filling digi histos for wheel and for RZ and RPhi
306  WheelHistos(id.wheel())->Fill((*digiIt).time(), theta, id.superlayer());
307  }
308 
309  } // for digis in layer
310 
311  meDoubleDigi_->Fill((float)num_digis_layer);
312 
313  } // for layers
314 
315  // cout << "num_digis " << num_digis << "mu digis " << num_mudigis << endl;
316 
317  if (num_musimhits != 0) {
318  meDigiEfficiencyMu_->Fill((float)num_mudigis / (float)num_musimhits);
319  meDigiEfficiency_->Fill((float)num_digis / (float)num_musimhits);
320  }
321 
322  meSimvsDigi_->Fill((float)num_musimhits, (float)num_digis);
323  // cout<<"--------------"<<endl;
324 }
325 
327  switch (abs(wheel)) {
328  case 0:
329  return hDigis_W0.get();
330 
331  case 1:
332  return hDigis_W1.get();
333 
334  case 2:
335  return hDigis_W2.get();
336 
337  default:
338  return nullptr;
339  }
340 }
344 
FastTimerService_cff.range
range
Definition: FastTimerService_cff.py:34
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
amptDefaultParameters_cff.mu
mu
Definition: amptDefaultParameters_cff.py:16
MuonDTDigis::DTWireIdMap
std::map< DTWireId, std::vector< const PSimHit * > > DTWireIdMap
Definition: MuonDTDigis.h:64
edm::Run
Definition: Run.h:45
MuonDTDigis.h
relativeConstraints.station
station
Definition: relativeConstraints.py:67
edm
HLT enums.
Definition: AlignableModifier.h:19
MuonDTDigis::analyze
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override
Definition: MuonDTDigis.cc:162
gather_cfg.cout
cout
Definition: gather_cfg.py:144
dqm::implementation::NavigatorBase::setCurrentFolder
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:32
FastTrackerRecHitCombiner_cfi.simHits
simHits
Definition: FastTrackerRecHitCombiner_cfi.py:5
DQMStore.h
edm::Handle< DTDigiCollection >
MuonDigiCollection::const_iterator
std::vector< DigiType >::const_iterator const_iterator
Definition: MuonDigiCollection.h:94
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
MakerMacros.h
edm::EventSetup::get
T get() const
Definition: EventSetup.h:87
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
DTWireId
Definition: DTWireId.h:12
edm::ESHandle< DTGeometry >
MuonDTDigis::~MuonDTDigis
~MuonDTDigis() override
Definition: MuonDTDigis.cc:37
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
Point3DBase< float, LocalTag >
DTLayerId
Definition: DTLayerId.h:12
edm::ParameterSet
Definition: ParameterSet.h:47
edm::get
T const & get(Event const &event, InputTag const &tag) noexcept(false)
Definition: Event.h:675
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
makeMuonMisalignmentScenario.wheel
wheel
Definition: makeMuonMisalignmentScenario.py:319
ModuleDef.h
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
edm::EventSetup
Definition: EventSetup.h:58
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
MuonDTDigis::MuonDTDigis
MuonDTDigis(const edm::ParameterSet &pset)
Definition: MuonDTDigis.cc:19
std
Definition: JetResolutionObject.h:76
dqm::implementation::IBooker::book2D
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:177
MuonDTDigis
Definition: MuonDTDigis.h:48
MuonDigiCollection::Range
std::pair< const_iterator, const_iterator > Range
Definition: MuonDigiCollection.h:95
dqm::implementation::IBooker
Definition: DQMStore.h:43
MuonDTDigis::WheelHistos
hDigis * WheelHistos(int wheel)
Definition: MuonDTDigis.cc:326
castor_dqm_sourceclient_file_cfg.path
path
Definition: castor_dqm_sourceclient_file_cfg.py:37
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
event
Definition: event.py:1
DigiContainerIterator
Definition: MuonDigiCollection.h:30
edm::Event
Definition: Event.h:73
MuonGeometryRecord
Definition: MuonGeometryRecord.h:34
edm::InputTag
Definition: InputTag.h:15
hit
Definition: SiStripHitEffFromCalibTree.cc:88
MuonDTDigis::bookHistograms
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: MuonDTDigis.cc:42
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
dqm::implementation::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98