CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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  // ----------------------
22  // Get the debug parameter for verbose output
23  verbose_ = pset.getUntrackedParameter<bool>("verbose",false);
24 
25  // the name of the Digi collection
26  SimHitToken_ = consumes<edm::PSimHitContainer>(pset.getParameter<edm::InputTag>("SimHitLabel"));
27 
28  // the name of the Digi collection
29  DigiToken_ = consumes<DTDigiCollection> (pset.getParameter<edm::InputTag>(("DigiLabel")));
30 
31  // ----------------------
32  // DQM ROOT output
33  outputFile_ = pset.getUntrackedParameter<std::string>("outputFile", "");
34  if ( outputFile_.size() != 0 ) {
35  LogInfo("OutputInfo") << " DT Muon Digis Task histograms will be saved to '"
36  << outputFile_.c_str() << "'";
37  } else {
38  LogInfo("OutputInfo") << " DT Muon Digis Task histograms will NOT be saved";
39  }
40 
41  string::size_type loc = outputFile_.find( ".root", 0 );
42  std::string outputFile_more_plots_;
43  if( loc != string::npos ) {
44  outputFile_more_plots_ = outputFile_.substr(0,loc)+"_more_plots.root";
45  } else {
46  outputFile_more_plots_ = " DTDigis_more_plots.root";
47  }
48 
49  // Please, uncomment next lines if you want a secong root file with additional histos
50 
51  // file_more_plots = new TFile(outputFile_more_plots_.c_str(),"RECREATE");
52  // file_more_plots->cd();
53  // if(file_more_plots->IsOpen()) cout<<"File for additional plots: " << outputFile_more_plots_ << " open!"<<endl;
54  // else cout<<"*** Error in opening file for additional plots ***"<<endl;
55 
56  hDigis_global = new hDigis("Global");
57  hDigis_W0 = new hDigis("Wheel0");
58  hDigis_W1 = new hDigis("Wheel1");
59  hDigis_W2 = new hDigis("Wheel2");
60  hAllHits = new hHits("AllHits");
61 
62  // End of comment . See more in Destructor (~MuonDTDigis) method
63 
64 
65 
66  // ----------------------
67  // get hold of back-end interface
68  dbe_ = 0;
70  if ( dbe_ ) {
71  if ( verbose_ ) {
72  dbe_->setVerbose(1);
73  } else {
74  dbe_->setVerbose(0);
75  }
76  }
77  if ( dbe_ ) {
78  if ( verbose_ ) dbe_->showDirStructure();
79  }
80 
81  // ----------------------
82 
83  meDigiTimeBox_ = 0;
84  meDigiTimeBox_wheel2m_ = 0;
85  meDigiTimeBox_wheel1m_ = 0;
86  meDigiTimeBox_wheel0_ = 0;
87  meDigiTimeBox_wheel1p_ = 0;
88  meDigiTimeBox_wheel2p_ = 0;
89  meDigiEfficiency_ = 0;
90  meDigiEfficiencyMu_ = 0;
91  meDoubleDigi_ = 0;
92  meSimvsDigi_ = 0;
93  meWire_DoubleDigi_ = 0;
94 
95  meMB1_sim_occup_ = 0;
96  meMB1_digi_occup_ = 0;
97  meMB2_sim_occup_ = 0;
98  meMB2_digi_occup_ = 0;
99  meMB3_sim_occup_ = 0;
100  meMB3_digi_occup_ = 0;
101  meMB4_sim_occup_ = 0;
102  meMB4_digi_occup_ = 0;
103 
104  //meDigiTimeBox_SL_ = 0;
105  meDigiHisto_ = 0;
106 
107 
108  // ----------------------
109  // We go
110  // ----------------------
111 
112  Char_t histo_n[100];
113  Char_t histo_t[100];
114 
115  if ( dbe_ ) {
116  dbe_->setCurrentFolder("MuonDTDigisV/DTDigiValidationTask");
117 
118  sprintf (histo_n, "DigiTimeBox" );
119  sprintf (histo_t, "Digi Time Box" );
120  meDigiTimeBox_ = dbe_->book1D(histo_n, histo_t, 1536,0,1200);
121 
122  sprintf (histo_n, "DigiTimeBox_wheel2m" );
123  sprintf (histo_t, "Digi Time Box wheel -2" );
124  meDigiTimeBox_wheel2m_ = dbe_->book1D(histo_n, histo_t, 384,0,1200);
125 
126  sprintf (histo_n, "DigiTimeBox_wheel1m" );
127  sprintf (histo_t, "Digi Time Box wheel -1" );
128  meDigiTimeBox_wheel1m_ = dbe_->book1D(histo_n, histo_t, 384,0,1200);
129 
130  sprintf (histo_n, "DigiTimeBox_wheel0" );
131  sprintf (histo_t, "Digi Time Box wheel 0" );
132  meDigiTimeBox_wheel0_ = dbe_->book1D(histo_n, histo_t, 384,0,1200);
133 
134  sprintf (histo_n, "DigiTimeBox_wheel1p" );
135  sprintf (histo_t, "Digi Time Box wheel 1" );
136  meDigiTimeBox_wheel1p_ = dbe_->book1D(histo_n, histo_t, 384,0,1200);
137 
138  sprintf (histo_n, "DigiTimeBox_wheel2p" );
139  sprintf (histo_t, "Digi Time Box wheel 2" );
140  meDigiTimeBox_wheel2p_ = dbe_->book1D(histo_n, histo_t, 384,0,1200);
141 
142  sprintf (histo_n, "DigiEfficiencyMu" );
143  sprintf (histo_t, "Ratio (#Digis Mu)/(#SimHits Mu)" );
144  meDigiEfficiencyMu_ = dbe_->book1D(histo_n, histo_t, 100, 0., 5.);
145 
146  sprintf (histo_n, "DigiEfficiency" );
147  sprintf (histo_t, "Ratio (#Digis)/(#SimHits)" );
148  meDigiEfficiency_ = dbe_->book1D(histo_n, histo_t, 100, 0., 5.);
149 
150  sprintf (histo_n, "Number_Digi_per_layer" );
151  sprintf (histo_t, "Number_Digi_per_layer" );
152  meDoubleDigi_ = dbe_->book1D(histo_n, histo_t, 10,0.,10.);
153 
154  sprintf (histo_n, "Number_simhit_vs_digi" );
155  sprintf (histo_t, "Number_simhit_vs_digi" );
156  meSimvsDigi_ = dbe_->book2D(histo_n, histo_t, 100, 0., 140., 100, 0., 140.);
157 
158  sprintf (histo_n, "Wire_Number_with_double_Digi" );
159  sprintf (histo_t, "Wire_Number_with_double_Digi" );
160  meWire_DoubleDigi_ = dbe_->book1D(histo_n, histo_t, 100,0.,100.);
161 
162  sprintf (histo_n, "Simhit_occupancy_MB1" );
163  sprintf (histo_t, "Simhit_occupancy_MB1" );
164  meMB1_sim_occup_ = dbe_->book1D(histo_n, histo_t, 55, 0., 55. );
165 
166  sprintf (histo_n, "Digi_occupancy_MB1" );
167  sprintf (histo_t, "Digi_occupancy_MB1" );
168  meMB1_digi_occup_ = dbe_->book1D(histo_n, histo_t, 55, 0., 55. );
169 
170  sprintf (histo_n, "Simhit_occupancy_MB2" );
171  sprintf (histo_t, "Simhit_occupancy_MB2" );
172  meMB2_sim_occup_ = dbe_->book1D(histo_n, histo_t, 63, 0., 63. );
173 
174  sprintf (histo_n, "Digi_occupancy_MB2" );
175  sprintf (histo_t, "Digi_occupancy_MB2" );
176  meMB2_digi_occup_ = dbe_->book1D(histo_n, histo_t, 63, 0., 63. );
177 
178  sprintf (histo_n, "Simhit_occupancy_MB3" );
179  sprintf (histo_t, "Simhit_occupancy_MB3" );
180  meMB3_sim_occup_ = dbe_->book1D(histo_n, histo_t, 75, 0., 75. );
181 
182  sprintf (histo_n, "Digi_occupancy_MB3" );
183  sprintf (histo_t, "Digi_occupancy_MB3" );
184  meMB3_digi_occup_ = dbe_->book1D(histo_n, histo_t, 75, 0., 75. );
185 
186  sprintf (histo_n, "Simhit_occupancy_MB4" );
187  sprintf (histo_t, "Simhit_occupancy_MB4" );
188  meMB4_sim_occup_ = dbe_->book1D(histo_n, histo_t, 99, 0., 99. );
189 
190  sprintf (histo_n, "Digi_occupancy_MB4" );
191  sprintf (histo_t, "Digi_occupancy_MB4" );
192  meMB4_digi_occup_ = dbe_->book1D(histo_n, histo_t, 99, 0., 99. );
193 
194  // sprintf (histo_n, "" );
195  // sprintf (histo_t, "" );
196 
197  /*
198  // Other option
199  string histoTag = "DigiTimeBox_slid_";
200  for ( int wheel = -2; wheel <= +2; ++wheel ) {
201  for ( int station = 1; station <= 4; ++station ) {
202  for ( int superLayer = 1; superLayer <= 3; ++superLayer ) {
203  string histoName = histoTag
204  + "_W" + wheel.str()
205  + "_St" + station.str()
206  + "_SL" + superLayer.str();
207  // the booking is not yet done
208  }
209  }
210  }
211  */
212  // Begona's option
213  char stringcham[40];
214  for ( int slnum = 1; slnum < 62; ++slnum ) {
215  sprintf(stringcham, "DigiTimeBox_slid_%d", slnum) ;
216  meDigiHisto_ = dbe_->book1D(stringcham, stringcham, 100,0,1200);
217  meDigiTimeBox_SL_.push_back(meDigiHisto_);
218  }
219  }
220 
221 }
222 
224 
225  // Uncomment these next lines if you want additional plots in a second .root file
226 
227  // file_more_plots->cd();
228 
229  // hDigis_global->Write();
230 
231  // hDigis_W0->Write();
232  // hDigis_W1->Write();
233  // hDigis_W2->Write();
234  // hAllHits->Write();
235 
236  // file_more_plots->Close();
237 
238  // End of comment.
239 
240  if ( outputFile_.size() != 0 && dbe_ ) dbe_->save(outputFile_);
241 
242  if(verbose_)
243  cout << "[MuonDTDigis] Destructor called" << endl;
244 
245 }
246 
247 //void MuonDTDigis::endJob(void){
248 
249 //}
250 
251 void MuonDTDigis::analyze(const Event & event, const EventSetup& eventSetup){
252 
253  if(verbose_)
254  cout << "--- [MuonDTDigis] Analysing Event: #Run: " << event.id().run()
255  << " #Event: " << event.id().event() << endl;
256 
257  // Get the DT Geometry
258  ESHandle<DTGeometry> muonGeom;
259  eventSetup.get<MuonGeometryRecord>().get(muonGeom);
260 
261  // Get the Digi collection from the event
262  Handle<DTDigiCollection> dtDigis;
263  event.getByToken(DigiToken_, dtDigis);
264 
265  // Get simhits
267  event.getByToken(SimHitToken_, simHits);
268 
269 
270  int num_mudigis;
271  int num_musimhits;
272  int num_digis;
273  int num_digis_layer;
274  int cham_num ;
275  int wire_touched;
276  num_digis = 0;
277  num_mudigis = 0;
278  num_musimhits = 0;
279  DTWireIdMap wireMap;
280 
281  for(vector<PSimHit>::const_iterator hit = simHits->begin();
282  hit != simHits->end(); hit++){
283  // Create the id of the wire, the simHits in the DT known also the wireId
284  DTWireId wireId(hit->detUnitId());
285  // cout << " PSimHits wire id " << wireId << " part type " << hit->particleType() << endl;
286 
287  // Fill the map
288  wireMap[wireId].push_back(&(*hit));
289 
290  LocalPoint entryP = hit->entryPoint();
291  LocalPoint exitP = hit->exitPoint();
292  int partType = hit->particleType();
293  if ( abs(partType) == 13 ) num_musimhits++;
294 
295  if ( wireId.station() == 1 && abs(partType) == 13 ) meMB1_sim_occup_->Fill(wireId.wire());
296  if ( wireId.station() == 2 && abs(partType) == 13 ) meMB2_sim_occup_->Fill(wireId.wire());
297  if ( wireId.station() == 3 && abs(partType) == 13 ) meMB3_sim_occup_->Fill(wireId.wire());
298  if ( wireId.station() == 4 && abs(partType) == 13 ) meMB4_sim_occup_->Fill(wireId.wire());
299 
300 
301  float path = (exitP-entryP).mag();
302  float path_x = fabs((exitP-entryP).x());
303 
304  hAllHits->Fill(entryP.x(),exitP.x(),
305  entryP.y(),exitP.y(),
306  entryP.z(),exitP.z(),
307  path , path_x,
308  partType, hit->processType(),
309  hit->pabs());
310 
311  }
312 
313  // cout << "num muon simhits " << num_musimhits << endl;
314 
316  for (detUnitIt=dtDigis->begin();
317  detUnitIt!=dtDigis->end();
318  ++detUnitIt){
319 
320  const DTLayerId& id = (*detUnitIt).first;
321  const DTDigiCollection::Range& range = (*detUnitIt).second;
322 
323  num_digis_layer = 0 ;
324  cham_num = 0 ;
325  wire_touched = 0;
326 
327  // Loop over the digis of this DetUnit
328  for (DTDigiCollection::const_iterator digiIt = range.first;
329  digiIt!=range.second;
330  ++digiIt){
331  // cout<<" Wire: "<<(*digiIt).wire()<<endl
332  // <<" digi time (ns): "<<(*digiIt).time()<<endl;
333 
334  num_digis++;
335  num_digis_layer++;
336  if (num_digis_layer > 1 )
337  {
338  if ( (*digiIt).wire() == wire_touched )
339  {
340  meWire_DoubleDigi_->Fill((*digiIt).wire());
341  // cout << "old & new wire " << wire_touched << " " << (*digiIt).wire() << endl;
342  }
343  }
344  wire_touched = (*digiIt).wire();
345 
346  meDigiTimeBox_->Fill((*digiIt).time());
347  if (id.wheel() == -2 ) meDigiTimeBox_wheel2m_->Fill((*digiIt).time());
348  if (id.wheel() == -1 ) meDigiTimeBox_wheel1m_->Fill((*digiIt).time());
349  if (id.wheel() == 0 ) meDigiTimeBox_wheel0_ ->Fill((*digiIt).time());
350  if (id.wheel() == 1 ) meDigiTimeBox_wheel1p_->Fill((*digiIt).time());
351  if (id.wheel() == 2 ) meDigiTimeBox_wheel2p_->Fill((*digiIt).time());
352 
353  // Superlayer number and fill histo with digi timebox
354 
355  cham_num = (id.wheel() +2)*12 + (id.station() -1)*3 + id.superlayer();
356  // cout << " Histo number " << cham_num << endl;
357 
358  meDigiTimeBox_SL_[cham_num]->Fill((*digiIt).time());
359 
360  // cout << " size de digis " << (*digiIt).size() << endl;
361 
362  DTWireId wireId(id,(*digiIt).wire());
363  if (wireId.station() == 1 ) meMB1_digi_occup_->Fill((*digiIt).wire());
364  if (wireId.station() == 2 ) meMB2_digi_occup_->Fill((*digiIt).wire());
365  if (wireId.station() == 3 ) meMB3_digi_occup_->Fill((*digiIt).wire());
366  if (wireId.station() == 4 ) meMB4_digi_occup_->Fill((*digiIt).wire());
367 
368  int mu=0;
369  float theta = 0;
370 
371  for(vector<const PSimHit*>::iterator hit = wireMap[wireId].begin();
372  hit != wireMap[wireId].end(); hit++)
373  if( abs((*hit)->particleType()) == 13){
374  theta = atan( (*hit)->momentumAtEntry().x()/ (-(*hit)->momentumAtEntry().z()) )*180/M_PI;
375  // cout<<"momentum x: "<<(*hit)->momentumAtEntry().x()<<endl
376  // <<"momentum z: "<<(*hit)->momentumAtEntry().z()<<endl
377  // <<"atan: "<<theta<<endl;
378  mu++;
379  }
380 
381  if( mu ) num_mudigis++;
382 
383  if(mu && theta){
384  hDigis_global->Fill((*digiIt).time(),theta,id.superlayer());
385  //filling digi histos for wheel and for RZ and RPhi
386  WheelHistos(id.wheel())->Fill((*digiIt).time(),theta,id.superlayer());
387  }
388 
389  }// for digis in layer
390 
391  meDoubleDigi_->Fill( (float)num_digis_layer );
392 
393  }// for layers
394 
395  //cout << "num_digis " << num_digis << "mu digis " << num_mudigis << endl;
396 
397  if (num_musimhits != 0) {
398  meDigiEfficiencyMu_->Fill( (float)num_mudigis/(float)num_musimhits );
399  meDigiEfficiency_->Fill( (float)num_digis/(float)num_musimhits );
400  }
401 
402  meSimvsDigi_->Fill( (float)num_musimhits, (float)num_digis ) ;
403  // cout<<"--------------"<<endl;
404 
405 }
406 
407 hDigis* MuonDTDigis::WheelHistos(int wheel){
408  switch(abs(wheel)){
409 
410  case 0: return hDigis_W0;
411 
412  case 1: return hDigis_W1;
413 
414  case 2: return hDigis_W2;
415 
416  default: return NULL;
417  }
418 }
422 
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
Definition: MuonDTDigis.cc:251
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:954
std::map< DTWireId, std::vector< const PSimHit * > > DTWireIdMap
Definition: MuonDTDigis.h:71
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Geom::Theta< T > theta() const
T y() const
Definition: PV3DBase.h:63
#define NULL
Definition: scimark2.h:8
DEFINE_FWK_MODULE(HiMixingModule)
uint16_t size_type
tuple path
else: Piece not in the list, fine.
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
T z() const
Definition: PV3DBase.h:64
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual ~MuonDTDigis()
Definition: MuonDTDigis.cc:223
const int mu
Definition: Constants.h:22
void setVerbose(unsigned level)
Definition: DQMStore.cc:631
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
#define M_PI
hDigis * WheelHistos(int wheel)
Definition: MuonDTDigis.cc:407
DQMStore * dbe_
const T & get() const
Definition: EventSetup.h:55
std::vector< DTDigi >::const_iterator const_iterator
tuple simHits
Definition: trackerHits.py:16
#define begin
Definition: vmac.h:30
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, const uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE", const bool resetMEsAfterWriting=false)
Definition: DQMStore.cc:2540
std::pair< const_iterator, const_iterator > Range
tuple cout
Definition: gather_cfg.py:121
void showDirStructure(void) const
Definition: DQMStore.cc:3332
Definition: DDAxes.h:10
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:1082
T x() const
Definition: PV3DBase.h:62
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:667
MuonDTDigis(const edm::ParameterSet &pset)
Definition: MuonDTDigis.cc:19