CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalRecHitsValidation.cc
Go to the documentation of this file.
4 
6  // DQM ROOT output
7  outputFile_ = conf.getUntrackedParameter<std::string>("outputFile", "myfile.root");
8 
9  if ( outputFile_.size() != 0 ) {
10  edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will be saved to '" << outputFile_.c_str() << "'";
11  } else {
12  edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will NOT be saved";
13  }
14 
15  nevtot = 0;
16 
17  dbe_ = 0;
18  // get hold of back-end interface
20 
21  Char_t histo[200];
22 
23  hcalselector_ = conf.getUntrackedParameter<std::string>("hcalselector", "all");
24  ecalselector_ = conf.getUntrackedParameter<std::string>("ecalselector", "yes");
25  eventype_ = conf.getUntrackedParameter<std::string>("eventype", "single");
26  sign_ = conf.getUntrackedParameter<std::string>("sign", "*");
27  mc_ = conf.getUntrackedParameter<std::string>("mc", "yes");
28  famos_ = conf.getUntrackedParameter<bool>("Famos", false);
29  useAllHistos_ = conf.getUntrackedParameter<bool>("useAllHistos", false);
30 
31  //Collections
32  tok_hbhe_ = consumes<HBHERecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HBHERecHitCollectionLabel"));
33  tok_hf_ = consumes<HFRecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HFRecHitCollectionLabel"));
34  tok_ho_ = consumes<HORecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HORecHitCollectionLabel"));
35 
36  // register for data access
37  tok_evt_ = consumes<edm::HepMCProduct>(edm::InputTag("generator"));
38  tok_EB_ = consumes<EBRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEB"));
39  tok_EE_ = consumes<EERecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
40  tok_hh_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits","HcalHits"));
41 
42  // std::cout << "*** famos_ = " << famos_ << std::endl;
43 
44  subdet_ = 5;
45  if (hcalselector_ == "noise") subdet_ = 0;
46  if (hcalselector_ == "HB" ) subdet_ = 1;
47  if (hcalselector_ == "HE" ) subdet_ = 2;
48  if (hcalselector_ == "HO" ) subdet_ = 3;
49  if (hcalselector_ == "HF" ) subdet_ = 4;
50  if (hcalselector_ == "all" ) subdet_ = 5;
51  if (hcalselector_ == "ZS" ) subdet_ = 6;
52 
53  etype_ = 1;
54  if (eventype_ == "multi") etype_ = 2;
55 
56  iz = 1;
57  if(sign_ == "-") iz = -1;
58  if(sign_ == "*") iz = 0;
59 
60  imc = 1;
61  if(mc_ == "no") imc = 0;
62 
63  if ( dbe_ ) {
64  dbe_->setCurrentFolder("HcalRecHitsV/HcalRecHitTask");
65 
66  // General counters (drawn)
67  sprintf (histo, "N_HB" );
68  Nhb = dbe_->book1D(histo, histo, 2600,0.,2600.);
69  sprintf (histo, "N_HE" );
70  Nhe = dbe_->book1D(histo, histo, 2600,0.,2600.);
71  sprintf (histo, "N_HO" );
72  Nho = dbe_->book1D(histo, histo, 2200,0.,2200.);
73  sprintf (histo, "N_HF" );
74  Nhf = dbe_->book1D(histo, histo, 1800,0., 1800.);
75 
76  // ZS
77  if(subdet_ == 6) {
78 
79  for (unsigned int i1 = 0; i1 < 82; i1++) {
80  for (unsigned int i2 = 0; i2 < 72; i2++) {
81  for (unsigned int i3 = 0; i3 < 4; i3++) {
82  for (unsigned int i4 = 0; i4 < 4; i4++) {
83  emap_min [i1][i2][i3][i4] = 99999.;
84  }
85  }
86  }
87  }
88 
89  //None of the ZS histos are drawn
90  if (useAllHistos_){
91  sprintf (histo, "ZSmin_map_depth1" );
92  map_depth1 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
93  sprintf (histo, "ZSmin_map_depth2" );
94  map_depth2 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
95  sprintf (histo, "ZSmin_map_depth3" );
96  map_depth3 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
97  sprintf (histo, "ZSmin_map_depth4" );
98  map_depth4 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
99 
100 
101  sprintf (histo, "ZS_Nreco_HB1" );
102  ZS_nHB1 = dbe_->book1D(histo, histo, 2500, 0., 2500.);
103  sprintf (histo, "ZS_Nreco_HB2" );
104  ZS_nHB2 = dbe_->book1D(histo, histo, 500, 0., 500.);
105  sprintf (histo, "ZS_Nreco_HE1" );
106  ZS_nHE1 = dbe_->book1D(histo, histo, 2000, 0., 2000.);
107  sprintf (histo, "ZS_Nreco_HE2" );
108  ZS_nHE2 = dbe_->book1D(histo, histo, 2000, 0., 2000.);
109  sprintf (histo, "ZS_Nreco_HE3" );
110  ZS_nHE3 = dbe_->book1D(histo, histo, 500, 0., 500.);
111  sprintf (histo, "ZS_Nreco_HO" );
112  ZS_nHO = dbe_->book1D(histo, histo, 2500, 0., 2500.);
113  sprintf (histo, "ZS_Nreco_HF1" );
114  ZS_nHF1 = dbe_->book1D(histo, histo, 1000, 0., 1000.);
115  sprintf (histo, "ZS_Nreco_HF2" );
116  ZS_nHF2 = dbe_->book1D(histo, histo, 1000, 0., 1000.);
117 
118  sprintf (histo, "ZSmin_simple1D_HB1" );
119  ZS_HB1 = dbe_->book1D(histo, histo,120, -2., 10.);
120  sprintf (histo, "ZSmin_simple1D_HB2" );
121  ZS_HB2 = dbe_->book1D(histo, histo,120, -2., 10.);
122  sprintf (histo, "ZSmin_simple1D_HE1" );
123  ZS_HE1 = dbe_->book1D(histo, histo,120, -2., 10.);
124  sprintf (histo, "ZSmin_simple1D_HE2" );
125  ZS_HE2 = dbe_->book1D(histo, histo,120, -2., 10.);
126  sprintf (histo, "ZSmin_simple1D_HE3" );
127  ZS_HE3 = dbe_->book1D(histo, histo,120, -2., 10.);
128  sprintf (histo, "ZSmin_simple1D_HO" );
129  ZS_HO = dbe_->book1D(histo, histo,120, -2., 10.);
130  sprintf (histo, "ZSmin_simple1D_HF1" );
131  ZS_HF1 = dbe_->book1D(histo, histo,200, -10., 10.);
132  sprintf (histo, "ZSmin_simple1D_HF2" );
133  ZS_HF2 = dbe_->book1D(histo, histo,200, -10., 10.);
134 
135  sprintf (histo, "ZSmin_sequential1D_HB1" );
136  ZS_seqHB1 = dbe_->book1D(histo, histo,2400, -1200., 1200.);
137  sprintf (histo, "ZSmin_sequential1D_HB2" );
138  ZS_seqHB2 = dbe_->book1D(histo, histo,2400, -1200., 1200.);
139  sprintf (histo, "ZSmin_sequential1D_HE1" );
140  ZS_seqHE1 = dbe_->book1D(histo, histo,4400, -2200., 2200.);
141  sprintf (histo, "ZSmin_sequential1D_HE2" );
142  ZS_seqHE2 = dbe_->book1D(histo, histo,4400, -2200., 2200.);
143  sprintf (histo, "ZSmin_sequential1D_HE3" );
144  ZS_seqHE3 = dbe_->book1D(histo, histo,4400, -2200., 2200.);
145  sprintf (histo, "ZSmin_sequential1D_HO" );
146  ZS_seqHO = dbe_->book1D(histo, histo,2400, -1200., 1200.);
147  sprintf (histo, "ZSmin_sequential1D_HF1" );
148  ZS_seqHF1 = dbe_->book1D(histo, histo,6000, -3000., 3000.);
149  sprintf (histo, "ZSmin_sequential1D_HF2" );
150  ZS_seqHF2 = dbe_->book1D(histo, histo,6000, -3000., 3000.);
151  }
152  }
153 
154  // ALL others, except ZS
155  else {
156 
157  sprintf (histo, "emap_depth1" );
158  emap_depth1 = dbe_->book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
159  sprintf (histo, "emap_depth2" );
160  emap_depth2 = dbe_->book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
161  sprintf (histo, "emap_depth3" );
162  emap_depth3 = dbe_->book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
163  sprintf (histo, "emap_depth4" );
164  emap_depth4 = dbe_->book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
165 
166  if (useAllHistos_){
167 
168  if (ecalselector_ == "yes") {
169  sprintf (histo, "map_ecal" );
170  map_ecal = dbe_->book2D(histo, histo, 70, -3.045, 3.045, 72, -3.1415926536, 3.1415926536);
171  }
172  }
173 
174  //The mean energy histos are drawn, but not the RMS or emean seq
175  sprintf (histo, "emean_vs_ieta_HB1" );
176  emean_vs_ieta_HB1 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s");
177  sprintf (histo, "emean_vs_ieta_HB2" );
178  emean_vs_ieta_HB2 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s");
179  sprintf (histo, "emean_vs_ieta_HE1" );
180  emean_vs_ieta_HE1 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10. ,2000., "s" );
181  sprintf (histo, "emean_vs_ieta_HE2" );
182  emean_vs_ieta_HE2 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s");
183  sprintf (histo, "emean_vs_ieta_HE3" );
184  emean_vs_ieta_HE3 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
185  sprintf (histo, "emean_vs_ieta_HO" );
186  emean_vs_ieta_HO = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
187  sprintf (histo, "emean_vs_ieta_HF1" );
188  emean_vs_ieta_HF1 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
189  sprintf (histo, "emean_vs_ieta_HF2" );
190  emean_vs_ieta_HF2 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
191 
192  if (useAllHistos_){
193  sprintf (histo, "RMS_vs_ieta_HB1" );
194  RMS_vs_ieta_HB1 = dbe_->book1D(histo, histo, 82, -41., 41.);
195  sprintf (histo, "RMS_vs_ieta_HB2" );
196  RMS_vs_ieta_HB2 = dbe_->book1D(histo, histo, 82, -41., 41.);
197  sprintf (histo, "RMS_vs_ieta_HE1" );
198  RMS_vs_ieta_HE1 = dbe_->book1D(histo, histo, 82, -41., 41.);
199  sprintf (histo, "RMS_vs_ieta_HE2" );
200  RMS_vs_ieta_HE2 = dbe_->book1D(histo, histo, 82, -41., 41.);
201  sprintf (histo, "RMS_vs_ieta_HE3" );
202  RMS_vs_ieta_HE3 = dbe_->book1D(histo, histo, 82, -41., 41.);
203  sprintf (histo, "RMS_vs_ieta_HO" );
204  RMS_vs_ieta_HO = dbe_->book1D(histo, histo, 82, -41., 41.);
205  sprintf (histo, "RMS_vs_ieta_HF1" );
206  RMS_vs_ieta_HF1 = dbe_->book1D(histo, histo, 82, -41., 41.);
207  sprintf (histo, "RMS_vs_ieta_HF2" );
208  RMS_vs_ieta_HF2 = dbe_->book1D(histo, histo, 82, -41., 41.);
209 
210  // Sequential emean and RMS
211  sprintf (histo, "emean_seq_HB1" );
212  emean_seqHB1 = dbe_->bookProfile(histo, histo, 2400, -1200., 1200., 2010, -10., 2000., "s" );
213  sprintf (histo, "emean_seq_HB2" );
214  emean_seqHB2 = dbe_->bookProfile(histo, histo, 2400, -1200., 1200., 2010, -10., 2000., "s" );
215  sprintf (histo, "emean_seq_HE1" );
216  emean_seqHE1 = dbe_->bookProfile(histo, histo, 4400, -2200., 2200., 2010, -10., 2000., "s" );
217  sprintf (histo, "emean_seq_HE2" );
218  emean_seqHE2 = dbe_->bookProfile(histo, histo, 4400, -2200., 2200., 2010, -10., 2000., "s" );
219  sprintf (histo, "emean_seq_HE3" );
220  emean_seqHE3 = dbe_->bookProfile(histo, histo, 4400, -2200., 2200., 2010, -10., 2000., "s" );
221  sprintf (histo, "emean_seq_HO" );
222  emean_seqHO = dbe_->bookProfile(histo, histo, 2400, -1200., 1200., 2010, -10., 2000., "s" );
223  sprintf (histo, "emean_seq_HF1" );
224  emean_seqHF1 = dbe_->bookProfile(histo, histo, 6000, -3000., 3000., 2010, -10., 2000., "s" );
225  sprintf (histo, "emean_seq_HF2" );
226  emean_seqHF2 = dbe_->bookProfile(histo, histo, 6000, -3000., 3000., 2010, -10., 2000., "s" );
227 
228  sprintf (histo, "RMS_seq_HB1" );
229  RMS_seq_HB1 = dbe_->book1D(histo, histo, 2400, -1200., 1200.);
230  sprintf (histo, "RMS_seq_HB2" );
231  RMS_seq_HB2 = dbe_->book1D(histo, histo, 2400, -1200., 1200.);
232  sprintf (histo, "RMS_seq_HE1" );
233  RMS_seq_HE1 = dbe_->book1D(histo, histo, 4400, -2200., 2200.);
234  sprintf (histo, "RMS_seq_HE2" );
235  RMS_seq_HE2 = dbe_->book1D(histo, histo, 4400, -2200., 2200.);
236  sprintf (histo, "RMS_seq_HE3" );
237  RMS_seq_HE3 = dbe_->book1D(histo, histo, 4400, -2200., 2200.);
238  sprintf (histo, "RMS_seq_HO" );
239  RMS_seq_HO = dbe_->book1D(histo, histo, 2400, -1200., 1200.);
240  sprintf (histo, "RMS_seq_HF1" );
241  RMS_seq_HF1 = dbe_->book1D(histo, histo, 6000, -3000., 3000.);
242  sprintf (histo, "RMS_seq_HF2" );
243  RMS_seq_HF2 = dbe_->book1D(histo, histo, 6000, -3000., 3000.);
244  }
245  // Occupancy
246  //The only occupancy histos drawn are occupancy vs. ieta
247  //but the maps are needed because this is where the latter are filled from
248  sprintf (histo, "occupancy_map_HB1" );
249  occupancy_map_HB1 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
250  sprintf (histo, "occupancy_map_HB2" );
251  occupancy_map_HB2 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
252  sprintf (histo, "occupancy_map_HE1" );
253  occupancy_map_HE1 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
254  sprintf (histo, "occupancy_map_HE2" );
255  occupancy_map_HE2 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
256  sprintf (histo, "occupancy_map_HE3" );
257  occupancy_map_HE3 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
258  sprintf (histo, "occupancy_map_HO" );
259  occupancy_map_HO = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
260  sprintf (histo, "occupancy_map_HF1" );
261  occupancy_map_HF1 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
262  sprintf (histo, "occupancy_map_HF2" );
263  occupancy_map_HF2 = dbe_->book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
264 
265  //These are drawn
266  sprintf (histo, "occupancy_vs_ieta_HB1" );
267  occupancy_vs_ieta_HB1 = dbe_->book1D(histo, histo, 82, -41., 41.);
268  sprintf (histo, "occupancy_vs_ieta_HB2" );
269  occupancy_vs_ieta_HB2 = dbe_->book1D(histo, histo, 82, -41., 41.);
270  sprintf (histo, "occupancy_vs_ieta_HE1" );
271  occupancy_vs_ieta_HE1 = dbe_->book1D(histo, histo, 82, -41., 41.);
272  sprintf (histo, "occupancy_vs_ieta_HE2" );
273  occupancy_vs_ieta_HE2 = dbe_->book1D(histo, histo, 82, -41., 41.);
274  sprintf (histo, "occupancy_vs_ieta_HE3" );
275  occupancy_vs_ieta_HE3 = dbe_->book1D(histo, histo, 82, -41., 41.);
276  sprintf (histo, "occupancy_vs_ieta_HO" );
277  occupancy_vs_ieta_HO = dbe_->book1D(histo, histo, 82, -41., 41.);
278  sprintf (histo, "occupancy_vs_ieta_HF1" );
279  occupancy_vs_ieta_HF1 = dbe_->book1D(histo, histo, 82, -41., 41.);
280  sprintf (histo, "occupancy_vs_ieta_HF2" );
281  occupancy_vs_ieta_HF2 = dbe_->book1D(histo, histo, 82, -41., 41.);
282 
283  //These are not
284  if (useAllHistos_){
285  sprintf (histo, "occ_sequential1D_HB1" );
286  occupancy_seqHB1 = dbe_->book1D(histo, histo,2400, -1200., 1200.);
287  sprintf (histo, "occ_sequential1D_HB2" );
288  occupancy_seqHB2 = dbe_->book1D(histo, histo,2400, -1200., 1200.);
289  sprintf (histo, "occ_sequential1D_HE1" );
290  occupancy_seqHE1 = dbe_->book1D(histo, histo,4400, -2200., 2200.);
291  sprintf (histo, "occ_sequential1D_HE2" );
292  occupancy_seqHE2 = dbe_->book1D(histo, histo,4400, -2200., 2200.);
293  sprintf (histo, "occ_sequential1D_HE3" );
294  occupancy_seqHE3 = dbe_->book1D(histo, histo,4400, -2200., 2200.);
295  sprintf (histo, "occ_sequential1D_HO" );
296  occupancy_seqHO = dbe_->book1D(histo, histo,2400, -1200., 1200.);
297  sprintf (histo, "occ_sequential1D_HF1" );
298  occupancy_seqHF1 = dbe_->book1D(histo, histo,6000, -3000., 3000.);
299  sprintf (histo, "occ_sequential1D_HF2" );
300  occupancy_seqHF2 = dbe_->book1D(histo, histo,6000, -3000., 3000.);
301  }
302 
303  //All status word histos except HF67 are drawn
304  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HB" ) ;
305  RecHit_StatusWord_HB = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
306 
307  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HE" ) ;
308  RecHit_StatusWord_HE = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
309 
310  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HF" ) ;
311  RecHit_StatusWord_HF = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
312 
313  if (useAllHistos_){
314  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HF67" ) ;
315  RecHit_StatusWord_HF67 = dbe_->book1D(histo, histo, 3 , 0.5, 3.5);
316  }
317  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HO" ) ;
318  RecHit_StatusWord_HO = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
319 
320  //Aux status word histos
321  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HB" ) ;
322  RecHit_Aux_StatusWord_HB = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
323 
324  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HE" ) ;
325  RecHit_Aux_StatusWord_HE = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
326 
327  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HF" ) ;
328  RecHit_Aux_StatusWord_HF = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
329 
330  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HO" ) ;
331  RecHit_Aux_StatusWord_HO = dbe_->book1D(histo, histo, 32 , -0.5, 31.5);
332 
333  //These are not drawn
334  if(imc !=0 && useAllHistos_) {
335  sprintf (histo, "map_econe_depth1" );
337  dbe_->book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
338  sprintf (histo, "map_econe_depth2" );
340  dbe_->book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
341  sprintf (histo, "map_econe_depth3" );
343  dbe_->book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
344  sprintf (histo, "map_econe_depth4" );
346  dbe_->book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
347  }
348  } // end-of (subdet_ =! 6)
349 
350  //======================= Now various cases one by one ===================
351 
352  //Histograms drawn for single pion scan
353  if(subdet_ != 0 && imc != 0) { // just not for noise
354  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths");
355  meEnConeEtaProfile = dbe_->bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
356 
357  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_E");
358  meEnConeEtaProfile_E = dbe_->bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
359 
360  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_EH");
361  meEnConeEtaProfile_EH = dbe_->bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
362  }
363  //The other cone profile, delta ieta/phi and noise histos are not drawn
364  if (useAllHistos_){
365  if(subdet_ != 0 && imc != 0) { // just not for noise
366 
367  // meEnConeEtaProfiel_depth1->Fill(eta_RecHit, HcalCone_d1);
368 
369  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth1");
370  meEnConeEtaProfile_depth1 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
371 
372  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth2");
373  meEnConeEtaProfile_depth2 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
374 
375  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth3");
376  meEnConeEtaProfile_depth3 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
377 
378  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth4");
379  meEnConeEtaProfile_depth4 = dbe_->bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
380 
381  }
382 
383  if(etype_ == 1 && subdet_ != 0) { // single part., not for noise
384 
385  sprintf (histo, "Delta_phi_cluster-MC");
386  meDeltaPhi = dbe_->book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
387 
388  sprintf (histo, "Delta_eta_cluster-MC");
389  meDeltaEta = dbe_->book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
390 
391  sprintf (histo, "Delta_phi_simcluster-MC");
392  meDeltaPhiS = dbe_->book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
393 
394  sprintf (histo, "Delta_eta_simcluster-MC");
395  meDeltaEtaS = dbe_->book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
396  }
397  // NOISE-specific
398 
399  if (hcalselector_ == "noise" ){
400 
401  sprintf (histo, "e_hb" ) ;
402  e_hb = dbe_->book1D(histo, histo,1000, -5., 5.);
403  sprintf (histo, "e_he" ) ;
404  e_he = dbe_->book1D(histo, histo,1000, -5., 5.);
405  sprintf (histo, "e_ho" ) ;
406  e_ho = dbe_->book1D(histo, histo,1000, -5., 5.);
407  sprintf (histo, "e_hfl" ) ;
408  e_hfl = dbe_->book1D(histo, histo,2000, -10., 10.);
409  sprintf (histo, "e_hfs" ) ;
410  e_hfs = dbe_->book1D(histo, histo,2000, -10., 10.);
411  }
412  }
413  // ************** HB **********************************
414  if (subdet_ == 1 || subdet_ == 5 ){
415 
416  //Only severity level, energy of rechits and overall HB timing histos are drawn
417  if (useAllHistos_){
418  if(etype_ == 1 && subdet_ == 1 ) {
419  if(imc != 0) {
420  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HB" ) ;
421  meNumRecHitsConeHB = dbe_->book1D(histo, histo, 100, 0., 100.);
422 
423  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HB" ) ;
424  meSumRecHitsEnergyConeHB = dbe_->book1D(histo,histo, 60 ,-20., 280.);
425  }
426 
427  sprintf (histo, "HcalRecHitTask_number_of_rechits_above_1GeV_HB");
428  meNumRecHitsThreshHB = dbe_->book1D(histo, histo, 30, 0., 30.);
429 
430  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HB" ) ;
431  meSumRecHitsEnergyHB = dbe_->book1D(histo,histo, 60 , -20., 280.);
432 
433  if (ecalselector_ == "yes") {
434  if(imc != 0) {
435  sprintf (histo, "HcalRecHitTask_number_of_ecalrechits_in_cone_HB");
436  meNumEcalRecHitsConeHB = dbe_->book1D(histo, histo, 300, 0., 300.);
437  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_in_cone_HB");
438  meEcalHcalEnergyConeHB = dbe_->book1D(histo,histo, 60 , -20., 280.);
439  }
440 
441  sprintf (histo, "HcalRecHitTask_energy_hcal_vs_ecal_HB");
442  meEnergyHcalVsEcalHB = dbe_->book2D(histo, histo, 300, 0., 150., 300, 0., 150.);
443  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_HB" ) ;
444  meEcalHcalEnergyHB = dbe_->book1D(histo,histo, 60 , -20., 280.);
445  }
446  }
447  }
448 
449  sprintf(histo, "HcalRecHitTask_severityLevel_HB");
450  sevLvl_HB = dbe_->book1D(histo, histo, 25, -0.5, 24.5);
451 
452  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HB" ) ;
453  meRecHitsEnergyHB = dbe_->book1D(histo, histo, 2010 , -10. , 2000.);
454 
455  sprintf (histo, "HcalRecHitTask_timing_HB" ) ;
456  meTimeHB = dbe_->book1D(histo, histo, 70, -48., 92.);
457 
458  //High, medium and low histograms to reduce RAM usage
459  sprintf (histo, "HcalRecHitTask_timing_vs_energy_Low_HB" ) ;
460  meTE_Low_HB = dbe_->book2D(histo, histo, 50, -5., 45., 70, -48., 92.);
461 
462  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HB" ) ;
463  meTE_HB = dbe_->book2D(histo, histo, 150, -5., 295., 70, -48., 92.);
464 
465  sprintf (histo, "HcalRecHitTask_timing_vs_energy_High_HB" ) ;
466  meTE_High_HB = dbe_->book2D(histo, histo, 150, -5., 2995., 70, -48., 92.);
467 
468  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HB" ) ;
469  meTEprofileHB_Low = dbe_->bookProfile(histo, histo, 50, -5., 45., 70, -48., 92.);
470 
471  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HB" ) ;
472  meTEprofileHB = dbe_->bookProfile(histo, histo, 150, -5., 295., 70, -48., 92.);
473 
474  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_High_HB" ) ;
475  meTEprofileHB_High = dbe_->bookProfile(histo, histo, 150, -5., 2995., 70, -48., 92.);
476 
477  //Timing by depth and rechits vs simhits are not drawn
478  if (useAllHistos_){
479  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HB_depth1" ) ;
480  meTE_HB1 = dbe_->book2D(histo, histo, 3000, -5., 2995., 70, -48., 92.);
481 
482  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HB_depth2" ) ;
483  meTE_HB2 = dbe_->book2D(histo, histo, 3000, -5., 2995., 70, -48., 92.);
484 
485  if(imc != 0) {
486  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HB");
487  meRecHitSimHitHB = dbe_->book2D(histo, histo, 120, 0., 1.2, 300, 0., 150.);
488  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HB");
489  meRecHitSimHitProfileHB = dbe_->bookProfile(histo, histo, 120, 0., 1.2, 500, 0., 500.);
490  }
491  }
492  }
493 
494  // ********************** HE ************************************
495  if ( subdet_ == 2 || subdet_ == 5 ){
496 
497  //None of these are drawn
498  if (useAllHistos_){
499  if(etype_ == 1 && subdet_ == 2 ) {
500 
501  if(imc != 0) {
502  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HE" ) ;
503  meNumRecHitsConeHE = dbe_->book1D(histo, histo, 100, 0., 100.);
504 
505  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HE" ) ;
506  meSumRecHitsEnergyConeHE = dbe_->book1D(histo,histo, 60 ,-20., 280.);
507  }
508 
509  sprintf (histo, "HcalRecHitTask_number_of_rechits_above_1GeV_HE");
510  meNumRecHitsThreshHE = dbe_->book1D(histo, histo, 30, 0., 30.);
511 
512  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HE" ) ;
513  meSumRecHitsEnergyHE = dbe_->book1D(histo,histo, 60 , -20., 280.);
514 
515  if (ecalselector_ == "yes") {
516  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_HE" ) ;
517  meEcalHcalEnergyHE = dbe_->book1D(histo,histo, 80, -20., 380.);
518 
519  sprintf (histo, "HcalRecHitTask_energy_hcal_vs_ecal_HE");
520  meEnergyHcalVsEcalHE = dbe_->book2D(histo, histo, 300, 0., 150., 300, 0., 150.);
521  if(imc != 0) {
522  sprintf (histo, "HcalRecHitTask_number_of_ecalrechits_in_cone_HE");
523  meNumEcalRecHitsConeHE = dbe_->book1D(histo, histo, 300, 0., 300.);
524  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_in_cone_HE");
525  meEcalHcalEnergyConeHE = dbe_->book1D(histo,histo, 60,-20., 280.);
526  }
527  }
528  }
529  }
530 
531  //Only severity level, energy of rechits and overall HB timing histos are drawn
532  sprintf(histo, "HcalRecHitTask_severityLevel_HE");
533  sevLvl_HE = dbe_->book1D(histo, histo, 25, -0.5, 24.5);
534 
535  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HE" ) ;
536  meRecHitsEnergyHE = dbe_->book1D(histo, histo, 510, -10., 500.);
537 
538  sprintf (histo, "HcalRecHitTask_timing_HE" ) ;
539  meTimeHE = dbe_->book1D(histo, histo, 70, -48., 92.);
540 
541  sprintf (histo, "HcalRecHitTask_timing_vs_energy_Low_HE" ) ;
542  meTE_Low_HE = dbe_->book2D(histo, histo, 80, -5., 75., 70, -48., 92.);
543 
544  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HE" ) ;
545  meTE_HE = dbe_->book2D(histo, histo, 200, -5., 395., 70, -48., 92.);
546 
547  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HE" ) ;
548  meTEprofileHE_Low = dbe_->bookProfile(histo, histo, 80, -5., 75., 70, -48., 92.);
549 
550  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HE" ) ;
551  meTEprofileHE = dbe_->bookProfile(histo, histo, 200, -5., 395., 70, -48., 92.);
552 
553  //Timing by depth and rechits vs simhits are not drawn
554  if (useAllHistos_){
555  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HE_depth1" ) ;
556  meTE_HE1 = dbe_->book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
557 
558  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HE_depth2" ) ;
559  meTE_HE2 = dbe_->book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
560 
561  if(imc != 0) {
562  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HE");
563  meRecHitSimHitHE = dbe_->book2D(histo, histo, 120, 0., 0.6, 300, 0., 150.);
564  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HE");
565  meRecHitSimHitProfileHE = dbe_->bookProfile(histo, histo, 120, 0., 0.6, 500, 0., 500.);
566  }
567  }
568 
569  }
570 
571  // ************** HO ****************************************
572  if ( subdet_ == 3 || subdet_ == 5 ){
573 
574  //Only severity level, energy of rechits and overall HB timing histos are drawn
575  if (useAllHistos_){
576  if(etype_ == 1 && subdet_ == 3) {
577  if (imc != 0) {
578  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HO" ) ;
579  meNumRecHitsConeHO = dbe_->book1D(histo, histo, 100, 0 , 100.);
580 
581  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HO" ) ;
582  meSumRecHitsEnergyConeHO = dbe_->book1D(histo,histo, 80 ,-20., 380.);
583  }
584 
585  sprintf (histo, "HcalRecHitTask_number_of_rechits_above_1GeV_HO");
586  meNumRecHitsThreshHO = dbe_->book1D(histo, histo, 100, 0., 100.);
587 
588  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HO" ) ;
589  meSumRecHitsEnergyHO = dbe_->book1D(histo,histo, 80 , -20., 380.);
590  }
591  }
592 
593  sprintf(histo, "HcalRecHitTask_severityLevel_HO");
594  sevLvl_HO = dbe_->book1D(histo, histo, 25, -0.5, 24.5);
595 
596  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HO" ) ;
597  meRecHitsEnergyHO = dbe_->book1D(histo, histo, 510 , -10. , 500.);
598 
599  sprintf (histo, "HcalRecHitTask_timing_HO" ) ;
600  meTimeHO = dbe_->book1D(histo, histo, 70, -48., 92.);
601 
602  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HO" ) ;
603  meTE_HO= dbe_->book2D(histo, histo, 60, -5., 55., 70, -48., 92.);
604 
605  sprintf (histo, "HcalRecHitTask_timing_vs_energy_High_HO" ) ;
606  meTE_High_HO= dbe_->book2D(histo, histo, 100, -5., 995., 70, -48., 92.);
607 
608  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HO" ) ;
609  meTEprofileHO = dbe_->bookProfile(histo, histo, 60, -5., 55., 70, -48., 92.);
610 
611  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_High_HO" ) ;
612  meTEprofileHO_High = dbe_->bookProfile(histo, histo, 100, -5., 995., 70, -48., 92.);
613 
614  //Rechits vs simhits are not drawn
615  if (useAllHistos_){
616  if(imc != 0) {
617  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HO");
618  meRecHitSimHitHO = dbe_->book2D(histo, histo, 150, 0., 1.5, 350, 0., 350.);
619  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HO");
620  meRecHitSimHitProfileHO = dbe_->bookProfile(histo, histo, 150, 0., 1.5, 500, 0., 500.);
621  }
622  }
623  }
624 
625  // ********************** HF ************************************
626  if ( subdet_ == 4 || subdet_ == 5 ){
627 
628  //Only severity level, energy of rechits and overall HB timing histos are drawn
629  if (useAllHistos_){
630  if(etype_ == 1 && subdet_ == 4) {
631 
632  if(imc != 0) {
633  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HF" ) ;
634  meNumRecHitsConeHF = dbe_->book1D(histo, histo, 30, 0 , 30.);
635 
636  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HF" ) ;
637  meSumRecHitsEnergyConeHF = dbe_->book1D(histo,histo,100, -20., 180.);
638 
639  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HFL" );
640  meSumRecHitsEnergyConeHFL = dbe_->book1D(histo,histo,100,-20., 180.);
641 
642  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HFS");
643  meSumRecHitsEnergyConeHFS = dbe_->book1D(histo,histo,100,-20., 180.);
644  }
645  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HF" ) ;
646  meSumRecHitsEnergyHF = dbe_->book1D(histo,histo, 80 , -20., 380.);
647  }
648  }
649 
650  sprintf(histo, "HcalRecHitTask_severityLevel_HF");
651  sevLvl_HF = dbe_->book1D(histo, histo, 25, -0.5, 24.5);
652 
653  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HF" ) ;
654  meRecHitsEnergyHF = dbe_->book1D(histo, histo, 1010 , -10. , 1000.);
655 
656  sprintf (histo, "HcalRecHitTask_timing_HF" ) ;
657  meTimeHF = dbe_->book1D(histo, histo, 70, -48., 92.);
658 
659  sprintf (histo, "HcalRecHitTask_timing_vs_energy_Low_HF" ) ;
660  meTE_Low_HF = dbe_->book2D(histo, histo, 100, -5., 195., 70, -48., 92.);
661 
662  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HF" ) ;
663  meTE_HF = dbe_->book2D(histo, histo, 200, -5., 995., 70, -48., 92.);
664 
665  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HF" ) ;
666  meTEprofileHF_Low = dbe_->bookProfile(histo, histo, 100, -5., 195., 70, -48., 92.);
667 
668  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HF" ) ;
669  meTEprofileHF = dbe_->bookProfile(histo, histo, 200, -5., 995., 70, -48., 92.);
670 
671  //Timing by L/S and rechits vs simhits are not drawn
672  if (useAllHistos_){
673  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HFL" ) ;
674  meTE_HFL = dbe_->book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
675 
676  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HFS" ) ;
677  meTE_HFS = dbe_->book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
678 
679  if(imc != 0) {
680  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HF");
681  meRecHitSimHitHF = dbe_->book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
682  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFL");
683  meRecHitSimHitHFL = dbe_->book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
684  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFS");
685  meRecHitSimHitHFS = dbe_->book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
686  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HF");
687  meRecHitSimHitProfileHF = dbe_->bookProfile(histo, histo, 50, 0., 50., 500, 0., 500.);
688  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFL");
689  meRecHitSimHitProfileHFL = dbe_->bookProfile(histo, histo, 50, 0., 50., 500, 0., 500.);
690  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFS");
691  meRecHitSimHitProfileHFS = dbe_->bookProfile(histo, histo, 50, 0., 50., 500, 0., 500.);
692  }
693  }
694  }
695  } //end-of if(_dbe)
696 
697 }
698 
699 
701 
703 
704 
706 
708 
709  using namespace edm;
710 
711  // cuts for each subdet_ector mimiking "Scheme B"
712  // double cutHB = 0.9, cutHE = 1.4, cutHO = 1.1, cutHFL = 1.2, cutHFS = 1.8;
713 
714  // energy in HCAL
715  double eHcal = 0.;
716  double eHcalCone = 0.;
717  double eHcalConeHB = 0.;
718  double eHcalConeHE = 0.;
719  double eHcalConeHO = 0.;
720  double eHcalConeHF = 0.;
721  double eHcalConeHFL = 0.;
722  double eHcalConeHFS = 0.;
723  // Total numbet of RecHits in HCAL, in the cone, above 1 GeV theshold
724  int nrechits = 0;
725  int nrechitsCone = 0;
726  int nrechitsThresh = 0;
727 
728  // energy in ECAL
729  double eEcal = 0.;
730  double eEcalB = 0.;
731  double eEcalE = 0.;
732  double eEcalCone = 0.;
733  int numrechitsEcal = 0;
734 
735  // MC info
736  double phi_MC = -999999.; // phi of initial particle from HepMC
737  double eta_MC = -999999.; // eta of initial particle from HepMC
738 
739  // HCAL energy around MC eta-phi at all depths;
740  double partR = 0.3;
741  double ehcal_coneMC_1 = 0.;
742  double ehcal_coneMC_2 = 0.;
743  double ehcal_coneMC_3 = 0.;
744  double ehcal_coneMC_4 = 0.;
745 
746  // Cone size for serach of the hottest HCAL cell around MC
747  double searchR = 1.0;
748  double eps = 0.001;
749 
750  // Single particle samples: actual eta-phi position of cluster around
751  // hottest cell
752  double etaHot = 99999.;
753  double phiHot = 99999.;
754 
755  // MC information
756 
757  // std::cout << "*** 1" << std::endl;
758 
759 
760  if(imc != 0) {
761 
763  ev.getByToken(tok_evt_,evtMC); // generator in late 310_preX
764  if (!evtMC.isValid()) {
765  std::cout << "no HepMCProduct found" << std::endl;
766  } else {
767  // std::cout << "*** source HepMCProduct found"<< std::endl;
768  }
769 
770  // MC particle with highest pt is taken as a direction reference
771  double maxPt = -99999.;
772  int npart = 0;
773  const HepMC::GenEvent * myGenEvent = evtMC->GetEvent();
774  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
775  p != myGenEvent->particles_end(); ++p ) {
776  double phip = (*p)->momentum().phi();
777  double etap = (*p)->momentum().eta();
778  // phi_MC = phip;
779  // eta_MC = etap;
780  double pt = (*p)->momentum().perp();
781  if(pt > maxPt) { npart++; maxPt = pt; phi_MC = phip; eta_MC = etap; }
782  }
783  // std::cout << "*** Max pT = " << maxPt << std::endl;
784 
785  }
786 
787  // std::cout << "*** 2" << std::endl;
788  // previously was: c.get<IdealGeometryRecord>().get (geometry);
789  c.get<CaloGeometryRecord>().get (geometry);
790 
791  // HCAL channel status map ****************************************
793  c.get<HcalChannelQualityRcd>().get( hcalChStatus );
794  theHcalChStatus = hcalChStatus.product();
795  // Assignment of severity levels **********************************
796  edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
797  c.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
798  theHcalSevLvlComputer = hcalSevLvlComputerHndl.product();
799 
800  // Fill working vectors of HCAL RecHits quantities (all of these are drawn)
801  fillRecHitsTmp(subdet_, ev);
802 
803  // HB
804  if( subdet_ ==5 || subdet_ == 1 ){
805  for(unsigned int iv=0; iv<hcalHBSevLvlVec.size(); iv++){
807  }
808  }
809  // HE
810  if( subdet_ ==5 || subdet_ == 2 ){
811  for(unsigned int iv=0; iv<hcalHESevLvlVec.size(); iv++){
813  }
814  }
815  // HO
816  if( subdet_ ==5 || subdet_ == 3 ){
817  for(unsigned int iv=0; iv<hcalHOSevLvlVec.size(); iv++){
819  }
820  }
821  // HF
822  if( subdet_ ==5 || subdet_ == 4 ){
823  for(unsigned int iv=0; iv<hcalHFSevLvlVec.size(); iv++){
825  }
826  }
827 
828  // std::cout << "*** 3" << std::endl;
829 
830 
831  //===========================================================================
832  // IN ALL other CASES : ieta-iphi maps
833  //===========================================================================
834 
835  // ECAL
836  if(ecalselector_ == "yes" && (subdet_ == 1 || subdet_ == 2 || subdet_ == 5)) {
838 
839 
840  ev.getByToken(tok_EB_, rhitEB);
841 
842  EcalRecHitCollection::const_iterator RecHit = rhitEB.product()->begin();
843  EcalRecHitCollection::const_iterator RecHitEnd = rhitEB.product()->end();
844 
845  for (; RecHit != RecHitEnd ; ++RecHit) {
846  EBDetId EBid = EBDetId(RecHit->id());
847 
848  const CaloCellGeometry* cellGeometry =
849  geometry->getSubdetectorGeometry (EBid)->getGeometry (EBid) ;
850  double eta = cellGeometry->getPosition ().eta () ;
851  double phi = cellGeometry->getPosition ().phi () ;
852  double en = RecHit->energy();
853  eEcal += en;
854  eEcalB += en;
855 
856  if (useAllHistos_) map_ecal->Fill(eta, phi, en);
857 
858  double r = dR(eta_MC, phi_MC, eta, phi);
859  if( r < partR) {
860  eEcalCone += en;
861  numrechitsEcal++;
862  }
863  }
864 
865 
867 
868  ev.getByToken(tok_EE_, rhitEE);
869 
870  RecHit = rhitEE.product()->begin();
871  RecHitEnd = rhitEE.product()->end();
872 
873  for (; RecHit != RecHitEnd ; ++RecHit) {
874  EEDetId EEid = EEDetId(RecHit->id());
875 
876  const CaloCellGeometry* cellGeometry =
877  geometry->getSubdetectorGeometry (EEid)->getGeometry (EEid) ;
878  double eta = cellGeometry->getPosition ().eta () ;
879  double phi = cellGeometry->getPosition ().phi () ;
880  double en = RecHit->energy();
881  eEcal += en;
882  eEcalE += en;
883 
884  if (useAllHistos_) map_ecal->Fill(eta, phi, en);
885 
886  double r = dR(eta_MC, phi_MC, eta, phi);
887  if( r < partR) {
888  eEcalCone += en;
889  numrechitsEcal++;
890  }
891  }
892  } // end of ECAL selection
893 
894 
895  // std::cout << "*** 4" << std::endl;
896 
897 
898  // Counting, including ZS items
899  // Filling HCAL maps ----------------------------------------------------
900  double maxE = -99999.;
901 
902  int nhb1 = 0;
903  int nhb2 = 0;
904  int nhe1 = 0;
905  int nhe2 = 0;
906  int nhe3 = 0;
907  int nho = 0;
908  int nhf1 = 0;
909  int nhf2 = 0;
910 
911  for (unsigned int i = 0; i < cen.size(); i++) {
912 
913  int sub = csub[i];
914  int depth = cdepth[i];
915  int ieta = cieta[i];
916  int iphi = ciphi[i];
917  double en = cen[i];
918  double eta = ceta[i];
919  double phi = cphi[i];
920  uint32_t stwd = cstwd[i];
921  uint32_t auxstwd = cauxstwd[i];
922  // double z = cz[i];
923 
924  int index = ieta * 72 + iphi; // for sequential histos
925 
926  /*
927  std::cout << "*** point 4-1" << " ieta, iphi, depth, sub = "
928  << ieta << ", " << iphi << ", " << depth << ", " << sub
929  << std::endl;
930  */
931 
932 
933  if( sub == 1 && depth == 1) nhb1++;
934  if( sub == 1 && depth == 2) nhb2++;
935  if( sub == 2 && depth == 1) nhe1++;
936  if( sub == 2 && depth == 2) nhe2++;
937  if( sub == 2 && depth == 3) nhe3++;
938  if( sub == 3 && depth == 4) nho++;
939  if( sub == 4 && depth == 1) nhf1++;
940  if( sub == 4 && depth == 2) nhf2++;
941 
942  if( subdet_ == 6) { // ZS specific
943  if( en < emap_min[ieta+41][iphi][depth-1][sub-1] )
944  emap_min[ieta+41][iphi][depth-1][sub-1] = en;
945  }
946 
947  double emin = 1.;
948  if(fabs(eta) > 3.) emin = 5.;
949 
950  double r = dR(eta_MC, phi_MC, eta, phi);
951  if( r < searchR ) { // search for hottest cell in a big cone
952  if(maxE < en && en > emin) {
953  maxE = en;
954  etaHot = eta;
955  phiHot = phi;
956  }
957  }
958 
959  /*
960  if(ieta == 27 ) {
961  std::cout << "*** ieta=28, iphi = " << iphi << " det = "
962  << sub << " depth = " << depth << std::endl;
963  }
964  */
965 
966  if( subdet_ != 6) {
967 
968  // std::cout << "*** 4-1" << std::endl;
969  //The emean_vs_ieta histos are drawn as well as the e_maps
970 
971 
972  // to distinguish HE and HF
973  if( depth == 1 || depth == 2 ) {
974  int ieta1 = ieta;
975  if(sub == 4) {
976  if (ieta1 < 0) ieta1--;
977  else ieta1++;
978  }
979  if (depth == 1) emap_depth1->Fill(double(ieta1), double(iphi), en);
980  if (depth == 2) emap_depth2->Fill(double(ieta1), double(iphi), en);
981  }
982 
983  if( depth == 3) emap_depth3->Fill(double(ieta), double(iphi), en);
984  if( depth == 4) emap_depth4->Fill(double(ieta), double(iphi), en);
985 
986  if (depth == 1 && sub == 1 ) {
987  emean_vs_ieta_HB1->Fill(double(ieta), en);
988  occupancy_map_HB1->Fill(double(ieta), double(iphi));
989  if(useAllHistos_){
990  emean_seqHB1->Fill(double(index), en);
991  }
992  }
993  if (depth == 2 && sub == 1) {
994  emean_vs_ieta_HB2->Fill(double(ieta), en);
995  occupancy_map_HB2->Fill(double(ieta), double(iphi));
996  if(useAllHistos_){
997  emean_seqHB2->Fill(double(index), en);
998  }
999  }
1000  if (depth == 1 && sub == 2) {
1001  emean_vs_ieta_HE1->Fill(double(ieta), en);
1002  occupancy_map_HE1->Fill(double(ieta), double(iphi));
1003  if(useAllHistos_){
1004  emean_seqHE1->Fill(double(index), en);
1005  }
1006  }
1007  if (depth == 2 && sub == 2) {
1008  emean_vs_ieta_HE2->Fill(double(ieta), en);
1009  occupancy_map_HE2->Fill(double(ieta), double(iphi));
1010  if(useAllHistos_){
1011  emean_seqHE2->Fill(double(index), en);
1012  }
1013  }
1014  if (depth == 3 && sub == 2) {
1015  emean_vs_ieta_HE3->Fill(double(ieta), en);
1016  occupancy_map_HE3->Fill(double(ieta), double(iphi));
1017  if(useAllHistos_){
1018  emean_seqHE3->Fill(double(index), en);
1019  }
1020  }
1021  if (depth == 4 ) {
1022  emean_vs_ieta_HO->Fill(double(ieta), en);
1023  occupancy_map_HO->Fill(double(ieta), double(iphi));
1024  if(useAllHistos_){
1025  emean_seqHO->Fill(double(index), en);
1026  }
1027  }
1028  if (depth == 1 && sub == 4) {
1029  emean_vs_ieta_HF1->Fill(double(ieta), en);
1030  occupancy_map_HF1->Fill(double(ieta), double(iphi));
1031  if(useAllHistos_){
1032  emean_seqHF1->Fill(double(index), en);
1033  }
1034  }
1035  if (depth == 2 && sub == 4) {
1036  emean_vs_ieta_HF2->Fill(double(ieta), en);
1037  occupancy_map_HF2->Fill(double(ieta), double(iphi));
1038  if(useAllHistos_){
1039  emean_seqHF2->Fill(double(index), en);
1040  }
1041  }
1042  }
1043 
1044 
1045  if( r < partR ) {
1046  if (depth == 1) ehcal_coneMC_1 += en;
1047  if (depth == 2) ehcal_coneMC_2 += en;
1048  if (depth == 3) ehcal_coneMC_3 += en;
1049  if (depth == 4) ehcal_coneMC_4 += en;
1050  }
1051 
1052  //32-bit status word
1053  uint32_t statadd;
1054  unsigned int isw67 = 0;
1055  for (unsigned int isw = 0; isw < 32; isw++){
1056  statadd = 0x1<<(isw);
1057  if (stwd & statadd){
1058  if (sub == 1) RecHit_StatusWord_HB->Fill(isw);
1059  else if (sub == 2) RecHit_StatusWord_HE->Fill(isw);
1060  else if (sub == 3) RecHit_StatusWord_HO->Fill(isw);
1061  else if (sub == 4){
1062  RecHit_StatusWord_HF->Fill(isw);
1063  if (isw == 6) isw67 += 1;
1064  if (isw == 7) isw67 += 2;
1065  }
1066  }
1067  }
1068  if (isw67 != 0 && useAllHistos_) RecHit_StatusWord_HF67->Fill(isw67); //This one is not drawn
1069 
1070  for (unsigned int isw =0; isw < 32; isw++){
1071  statadd = 0x1<<(isw);
1072  if( auxstwd & statadd ){
1073  if (sub == 1) RecHit_Aux_StatusWord_HB->Fill(isw);
1074  else if (sub == 2) RecHit_Aux_StatusWord_HE->Fill(isw);
1075  else if (sub == 3) RecHit_Aux_StatusWord_HO->Fill(isw);
1076  else if (sub == 4) RecHit_Aux_StatusWord_HF->Fill(isw);
1077  }
1078 
1079  }
1080 
1081  }
1082 
1083  // std::cout << "*** 4-2" << std::endl;
1084 
1085  if( subdet_ == 6 && useAllHistos_) { // ZS plots; not drawn
1086  ZS_nHB1->Fill(double(nhb1));
1087  ZS_nHB2->Fill(double(nhb2));
1088  ZS_nHE1->Fill(double(nhe1));
1089  ZS_nHE2->Fill(double(nhe2));
1090  ZS_nHE3->Fill(double(nhe3));
1091  ZS_nHO ->Fill(double(nho));
1092  ZS_nHF1->Fill(double(nhf1));
1093  ZS_nHF2->Fill(double(nhf2));
1094  }
1095  else{
1096  Nhb->Fill(double(nhb1 + nhb2));
1097  Nhe->Fill(double(nhe1 + nhe2 + nhe3));
1098  Nho->Fill(double(nho));
1099  Nhf->Fill(double(nhf1 + nhf2));
1100 
1101  //These are not drawn
1102  if(imc != 0 && useAllHistos_) {
1103  map_econe_depth1->Fill(eta_MC, phi_MC, ehcal_coneMC_1);
1104  map_econe_depth2->Fill(eta_MC, phi_MC, ehcal_coneMC_2);
1105  map_econe_depth3->Fill(eta_MC, phi_MC, ehcal_coneMC_3);
1106  map_econe_depth4->Fill(eta_MC, phi_MC, ehcal_coneMC_4);
1107  }
1108  }
1109 
1110  // std::cout << "*** 5" << std::endl;
1111 
1112 
1113  // NOISE =================================================================
1114  //Not drawn
1115  if (hcalselector_ == "noise" && useAllHistos_) {
1116  for (unsigned int i = 0; i < cen.size(); i++) {
1117 
1118  int sub = csub[i];
1119  int depth = cdepth[i];
1120  double en = cen[i];
1121 
1122  if (sub == 1) e_hb->Fill(en);
1123  if (sub == 2) e_he->Fill(en);
1124  if (sub == 3) e_ho->Fill(en);
1125  if (sub == 4) {
1126  if(depth == 1)
1127  e_hfl->Fill(en);
1128  else
1129  e_hfs->Fill(en);
1130  }
1131  }
1132  }
1133 
1134  //===========================================================================
1135  // SUBSYSTEMS,
1136  //===========================================================================
1137 
1138  else if ((subdet_ != 6) && (subdet_ != 0)) {
1139 
1140  // std::cout << "*** 6" << std::endl;
1141 
1142 
1143  double clusEta = 999.;
1144  double clusPhi = 999.;
1145  double clusEn = 0.;
1146 
1147  double HcalCone_d1 = 0.;
1148  double HcalCone_d2 = 0.;
1149  double HcalCone_d3 = 0.;
1150  double HcalCone_d4 = 0.;
1151  double HcalCone = 0.;
1152 
1153  int ietaMax1 = 9999;
1154  int ietaMax2 = 9999;
1155  int ietaMax3 = 9999;
1156  int ietaMax4 = 9999;
1157  int ietaMax = 9999;
1158  double enMax1 = -9999.;
1159  double enMax2 = -9999.;
1160  double enMax3 = -9999.;
1161  double enMax4 = -9999.;
1162  // double enMax = -9999.;
1163  double etaMax = 9999.;
1164 
1165  /*
1166  std::cout << "*** point 5-1" << " eta_MC, phi_MC etaHot, phiHot = "
1167  << eta_MC << ", " << phi_MC << " "
1168  << etaHot << ", " << phiHot
1169  << std::endl;
1170  */
1171 
1172  // CYCLE over cells ====================================================
1173 
1174  for (unsigned int i = 0; i < cen.size(); i++) {
1175  int sub = csub[i];
1176  int depth = cdepth[i];
1177  double eta = ceta[i];
1178  double phi = cphi[i];
1179  double en = cen[i];
1180  double t = ctime[i];
1181  int ieta = cieta[i];
1182 
1183  double rhot = dR(etaHot, phiHot, eta, phi);
1184  if(rhot < partR && en > 1.) {
1185  clusEta = (clusEta * clusEn + eta * en)/(clusEn + en);
1186  clusPhi = phi12(clusPhi, clusEn, phi, en);
1187  clusEn += en;
1188  }
1189 
1190  nrechits++;
1191  eHcal += en;
1192  if(en > 1. ) nrechitsThresh++;
1193 
1194  double r = dR(eta_MC, phi_MC, eta, phi);
1195  if( r < partR ){
1196  if(sub == 1) eHcalConeHB += en;
1197  if(sub == 2) eHcalConeHE += en;
1198  if(sub == 3) eHcalConeHO += en;
1199  if(sub == 4) {
1200  eHcalConeHF += en;
1201  if (depth == 1) eHcalConeHFL += en;
1202  else eHcalConeHFS += en;
1203  }
1204  eHcalCone += en;
1205  nrechitsCone++;
1206 
1207  // search for most energetic cell at the given depth in the cone
1208  if(depth == 1) {
1209  HcalCone_d1 += en;
1210  if(enMax1 < en) {
1211  enMax1 = en;
1212  ietaMax1 = ieta;
1213  }
1214  }
1215  if(depth == 2) {
1216  HcalCone_d2 += en;
1217  if(enMax2 < en) {
1218  enMax2 = en;
1219  ietaMax2 = ieta;
1220  }
1221  }
1222  if(depth == 3) {
1223  HcalCone_d3 += en;
1224  if(enMax3 < en) {
1225  enMax3 = en;
1226  ietaMax3 = ieta;
1227  }
1228  }
1229  if(depth == 4) {
1230  HcalCone_d4 += en;
1231  if(enMax4 < en) {
1232  enMax4 = en;
1233  ietaMax4 = ieta;
1234  }
1235  }
1236 
1237  if(depth != 4) {
1238  HcalCone += en;
1239  }
1240 
1241 
1242  // regardless of the depths (but excluding HO), just hottest cell
1243  /*
1244  if(depth != 4) {
1245  if(enMax < en) {
1246  enMax = en;
1247  ietaMax = ieta;
1248  }
1249  }
1250  */
1251 
1252  // alternative: ietamax -> closest to MC eta !!!
1253  float eta_diff = fabs(eta_MC - eta);
1254  if(eta_diff < etaMax) {
1255  etaMax = eta_diff;
1256  ietaMax = ieta;
1257  }
1258  }
1259 
1260  //The energy and overall timing histos are drawn while
1261  //the ones split by depth are not
1262  if(sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
1263  meTimeHB->Fill(t);
1264  meRecHitsEnergyHB->Fill(en);
1265 
1266  meTE_Low_HB->Fill( en, t);
1267  meTE_HB->Fill( en, t);
1268  meTE_High_HB->Fill( en, t);
1269  meTEprofileHB_Low->Fill(en, t);
1270  meTEprofileHB->Fill(en, t);
1271  meTEprofileHB_High->Fill(en, t);
1272 
1273  if (useAllHistos_){
1274  if (depth == 1) meTE_HB1->Fill( en, t);
1275  else if (depth == 2) meTE_HB2->Fill( en, t);
1276  }
1277  }
1278  if(sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
1279  meTimeHE->Fill(t);
1280  meRecHitsEnergyHE->Fill(en);
1281 
1282  meTE_Low_HE->Fill( en, t);
1283  meTE_HE->Fill( en, t);
1284  meTEprofileHE_Low->Fill(en, t);
1285  meTEprofileHE->Fill(en, t);
1286 
1287  if (useAllHistos_){
1288  if (depth == 1) meTE_HE1->Fill( en, t);
1289  else if (depth == 2) meTE_HE2->Fill( en, t);
1290  }
1291  }
1292  if(sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
1293  meTimeHF->Fill(t);
1294  meRecHitsEnergyHF->Fill(en);
1295 
1296  meTE_Low_HF->Fill(en, t);
1297  meTE_HF->Fill(en, t);
1298  meTEprofileHF_Low->Fill(en, t);
1299  meTEprofileHF->Fill(en, t);
1300 
1301  if (useAllHistos_){
1302  if (depth == 1) meTE_HFL->Fill( en, t);
1303  else meTE_HFS->Fill( en, t);
1304  }
1305  }
1306  if(sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
1307  meTimeHO->Fill(t);
1308  meRecHitsEnergyHO->Fill(en);
1309 
1310  meTE_HO->Fill( en, t);
1311  meTE_High_HO->Fill( en, t);
1312  meTEprofileHO->Fill(en, t);
1313  meTEprofileHO_High->Fill(en, t);
1314  }
1315  }
1316 
1317  if(imc != 0) {
1318  //Cone by depth are not drawn, the others are used for pion scan
1319  if (useAllHistos_){
1320  meEnConeEtaProfile_depth1->Fill(double(ietaMax1), HcalCone_d1);
1321  meEnConeEtaProfile_depth2->Fill(double(ietaMax2), HcalCone_d2);
1322  meEnConeEtaProfile_depth3->Fill(double(ietaMax3), HcalCone_d3);
1323  meEnConeEtaProfile_depth4->Fill(double(ietaMax4), HcalCone_d4);
1324  }
1325  meEnConeEtaProfile ->Fill(double(ietaMax), HcalCone); //
1326  meEnConeEtaProfile_E ->Fill(double(ietaMax), eEcalCone);
1327  meEnConeEtaProfile_EH ->Fill(double(ietaMax), HcalCone+eEcalCone);
1328  }
1329 
1330  // std::cout << "*** 7" << std::endl;
1331 
1332 
1333  // Single particle samples ONLY ! ======================================
1334  // Fill up some histos for "integrated" subsustems.
1335  // These are not drawn
1336  if(etype_ == 1 && useAllHistos_) {
1337 
1338  /*
1339  std::cout << "*** point 7-1" << " eta_MC, phi_MC clusEta, clusPhi = "
1340  << eta_MC << ", " << phi_MC << " "
1341  << clusEta << ", " << clusPhi
1342  << std::endl;
1343  */
1344 
1345  double phidev = dPhiWsign(clusPhi, phi_MC);
1346  meDeltaPhi->Fill(eta_MC, phidev);
1347  double etadev = clusEta - eta_MC;
1348  meDeltaEta->Fill(eta_MC, etadev);
1349 
1350  if(subdet_ == 1) {
1351  meSumRecHitsEnergyHB->Fill(eHcal);
1352  if(imc != 0) meSumRecHitsEnergyConeHB->Fill(eHcalConeHB);
1353  if(imc != 0) meNumRecHitsConeHB->Fill(double(nrechitsCone));
1354  meNumRecHitsThreshHB->Fill(double(nrechitsThresh));
1355  }
1356 
1357  if(subdet_ == 2) {
1358  meSumRecHitsEnergyHE->Fill(eHcal);
1359  if(imc != 0) meSumRecHitsEnergyConeHE->Fill(eHcalConeHE);
1360  if(imc != 0) meNumRecHitsConeHE->Fill(double(nrechitsCone));
1361  meNumRecHitsThreshHE->Fill(double(nrechitsThresh));
1362  }
1363 
1364  if(subdet_ == 3) {
1365  meSumRecHitsEnergyHO->Fill(eHcal);
1366  if(imc != 0) meSumRecHitsEnergyConeHO->Fill(eHcalConeHO);
1367  if(imc != 0) meNumRecHitsConeHO->Fill(double(nrechitsCone));
1368  meNumRecHitsThreshHO->Fill(double(nrechitsThresh));
1369  }
1370 
1371  if(subdet_ == 4) {
1372  if(eHcalConeHF > eps ) {
1373  meSumRecHitsEnergyHF ->Fill(eHcal);
1374  if(imc != 0) {
1375  meSumRecHitsEnergyConeHF ->Fill(eHcalConeHF);
1376  meNumRecHitsConeHF->Fill(double(nrechitsCone));
1377  meSumRecHitsEnergyConeHFL ->Fill(eHcalConeHFL);
1378  meSumRecHitsEnergyConeHFS ->Fill(eHcalConeHFS);
1379  }
1380  }
1381  }
1382 
1383  // std::cout << "*** 8" << std::endl;
1384 
1385 
1386  // Also combine with ECAL if needed
1387  if(subdet_ == 1 && ecalselector_ == "yes") {
1388 
1389  /*
1390  std::cout << "*** point 8-1"
1391  << " eEcalB " << eEcalB << " eHcal " << eHcal
1392  << " eEcalCone " << eEcalCone << " eHcalCone "
1393  << eHcalCone
1394  << " numrechitsEcal " << numrechitsEcal
1395  << std::endl;
1396 
1397  */
1398 
1399  meEcalHcalEnergyHB->Fill(eEcalB+eHcal);
1400  meEcalHcalEnergyConeHB->Fill(eEcalCone+eHcalCone);
1401  meNumEcalRecHitsConeHB->Fill(double(numrechitsEcal));
1402 
1403  }
1404 
1405  if(subdet_ == 2 && ecalselector_ == "yes"){
1406 
1407  /*
1408  std::cout << "*** point 8-2a"
1409  << " eEcalE " << eEcalE << " eHcal " << eHcal
1410  << " eEcalCone " << eEcalCone << " eHcalCone "
1411  << eHcalCone
1412  << " numrechitsEcal " << numrechitsEcal
1413  << std::endl;
1414  */
1415 
1416  meEcalHcalEnergyHE->Fill(eEcalE+eHcal);
1417  if(imc != 0) meEcalHcalEnergyConeHE->Fill(eEcalCone+eHcalCone);
1418  if(imc != 0) meNumEcalRecHitsConeHE->Fill(double(numrechitsEcal));
1419  }
1420 
1421  // Banana plots finally
1422  if(imc != 0) {
1423  if(subdet_ == 1 && ecalselector_ == "yes")
1424  meEnergyHcalVsEcalHB -> Fill(eEcalCone,eHcalCone);
1425  if(subdet_ == 2 && ecalselector_ == "yes")
1426  meEnergyHcalVsEcalHE -> Fill(eEcalCone,eHcalCone);
1427  }
1428  }
1429  }
1430  // std::cout << "*** 9" << std::endl;
1431 
1432 
1433  //===========================================================================
1434  // Getting SimHits
1435  //===========================================================================
1436 
1437  if(subdet_ > 0 && subdet_ < 6 && imc !=0 && !famos_ ) { // not noise
1438 
1439  double maxES = -9999.;
1440  double etaHotS = 1000.;
1441  double phiHotS = 1000.;
1442 
1444  ev.getByToken(tok_hh_,hcalHits);
1445  const PCaloHitContainer * SimHitResult = hcalHits.product () ;
1446 
1447  double enSimHits = 0.;
1448  double enSimHitsHB = 0.;
1449  double enSimHitsHE = 0.;
1450  double enSimHitsHO = 0.;
1451  double enSimHitsHF = 0.;
1452  double enSimHitsHFL = 0.;
1453  double enSimHitsHFS = 0.;
1454  // sum of SimHits in the cone
1455 
1456  for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin () ; SimHits != SimHitResult->end(); ++SimHits) {
1457  HcalDetId cell(SimHits->id());
1458  int sub = cell.subdet();
1459  const CaloCellGeometry* cellGeometry =
1460  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1461  double etaS = cellGeometry->getPosition().eta () ;
1462  double phiS = cellGeometry->getPosition().phi () ;
1463  double en = SimHits->energy();
1464 
1465  double emin = 0.01;
1466  if(fabs(etaS) > 3.) emin = 1.;
1467 
1468  double r = dR(eta_MC, phi_MC, etaS, phiS);
1469  if( r < searchR ) { // search for hottest cell in a big cone
1470  if(maxES < en && en > emin ) {
1471  maxES = en;
1472  etaHotS = etaS;
1473  phiHotS = phiS;
1474  }
1475  }
1476 
1477  if ( r < partR ){ // just energy in the small cone
1478  enSimHits += en;
1479  if(sub == 1) enSimHitsHB += en;
1480  if(sub == 2) enSimHitsHE += en;
1481  if(sub == 3) enSimHitsHO += en;
1482  if(sub == 4) {
1483  enSimHitsHF += en;
1484  int depth = cell.depth();
1485  if(depth == 1) enSimHitsHFL += en;
1486  else enSimHitsHFS += en;
1487  }
1488  }
1489  }
1490 
1491 
1492  // Second look over SimHits: cluster finding
1493 
1494  double clusEta = 999.;
1495  double clusPhi = 999.;
1496  double clusEn = 0.;
1497 
1498  for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin () ; SimHits != SimHitResult->end(); ++SimHits) {
1499  HcalDetId cell(SimHits->id());
1500 
1501  const CaloCellGeometry* cellGeometry =
1502  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1503  double etaS = cellGeometry->getPosition().eta () ;
1504  double phiS = cellGeometry->getPosition().phi () ;
1505  double en = SimHits->energy();
1506 
1507  double emin = 0.01;
1508  if(fabs(etaS) > 3.) emin = 1.;
1509 
1510  double rhot = dR(etaHotS, phiHotS, etaS, phiS);
1511  if(rhot < partR && en > emin) {
1512  clusEta = (clusEta * clusEn + etaS * en)/(clusEn + en);
1513  clusPhi = phi12(clusPhi, clusEn, phiS, en);
1514  clusEn += en;
1515  }
1516  }
1517 
1518  // SimHits cluster deviation from MC (eta, phi)
1519  // These are not drawn
1520  if (useAllHistos_){
1521  if(etype_ == 1) {
1522  double phidev = dPhiWsign(clusPhi, phi_MC);
1523  meDeltaPhiS->Fill(eta_MC, phidev);
1524  double etadev = clusEta - eta_MC;
1525  meDeltaEtaS->Fill(eta_MC, etadev);
1526  }
1527  // Now some histos with SimHits
1528 
1529  if(subdet_ == 4 || subdet_ == 5) {
1530  if(eHcalConeHF > eps) {
1531  meRecHitSimHitHF->Fill( enSimHitsHF, eHcalConeHF );
1532  meRecHitSimHitProfileHF->Fill( enSimHitsHF, eHcalConeHF);
1533 
1534  meRecHitSimHitHFL->Fill( enSimHitsHFL, eHcalConeHFL );
1535  meRecHitSimHitProfileHFL->Fill( enSimHitsHFL, eHcalConeHFL);
1536  meRecHitSimHitHFS->Fill( enSimHitsHFS, eHcalConeHFS );
1537  meRecHitSimHitProfileHFS->Fill( enSimHitsHFS, eHcalConeHFS);
1538  }
1539  }
1540  if(subdet_ == 1 || subdet_ == 5) {
1541  meRecHitSimHitHB->Fill( enSimHitsHB,eHcalConeHB );
1542  meRecHitSimHitProfileHB->Fill( enSimHitsHB,eHcalConeHB);
1543  }
1544  if(subdet_ == 2 || subdet_ == 5) {
1545  meRecHitSimHitHE->Fill( enSimHitsHE,eHcalConeHE );
1546  meRecHitSimHitProfileHE->Fill( enSimHitsHE,eHcalConeHE);
1547  }
1548  if(subdet_ == 3 || subdet_ == 5) {
1549  meRecHitSimHitHO->Fill( enSimHitsHO,eHcalConeHO );
1550  meRecHitSimHitProfileHO->Fill( enSimHitsHO,eHcalConeHO);
1551  }
1552  }
1553  }
1554 
1555  nevtot++;
1556 }
1557 
1558 
1561 
1562  using namespace edm;
1563 
1564 
1565  // initialize data vectors
1566  csub.clear();
1567  cen.clear();
1568  ceta.clear();
1569  cphi.clear();
1570  ctime.clear();
1571  cieta.clear();
1572  ciphi.clear();
1573  cdepth.clear();
1574  cz.clear();
1575  cstwd.clear();
1576  cauxstwd.clear();
1577  hcalHBSevLvlVec.clear();
1578  hcalHESevLvlVec.clear();
1579  hcalHFSevLvlVec.clear();
1580  hcalHOSevLvlVec.clear();
1581 
1582  if( subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1583 
1584  //HBHE
1586  ev.getByToken(tok_hbhe_, hbhecoll);
1587 
1588  for (HBHERecHitCollection::const_iterator j=hbhecoll->begin(); j != hbhecoll->end(); j++) {
1589  HcalDetId cell(j->id());
1590  const CaloCellGeometry* cellGeometry =
1591  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1592  double eta = cellGeometry->getPosition().eta () ;
1593  double phi = cellGeometry->getPosition().phi () ;
1594  double zc = cellGeometry->getPosition().z ();
1595  int sub = cell.subdet();
1596  int depth = cell.depth();
1597  int inteta = cell.ieta();
1598  if(inteta > 0) inteta -= 1;
1599  int intphi = cell.iphi()-1;
1600  double en = j->energy();
1601  double t = j->time();
1602  int stwd = j->flags();
1603  int auxstwd = j->aux();
1604 
1605  int serivityLevel = hcalSevLvl( (CaloRecHit*) &*j );
1606  if( cell.subdet()==HcalBarrel ){
1607  hcalHBSevLvlVec.push_back(serivityLevel);
1608  }else if (cell.subdet()==HcalEndcap ){
1609  hcalHESevLvlVec.push_back(serivityLevel);
1610  }
1611 
1612  if((iz > 0 && eta > 0.) || (iz < 0 && eta <0.) || iz == 0) {
1613 
1614  csub.push_back(sub);
1615  cen.push_back(en);
1616  ceta.push_back(eta);
1617  cphi.push_back(phi);
1618  ctime.push_back(t);
1619  cieta.push_back(inteta);
1620  ciphi.push_back(intphi);
1621  cdepth.push_back(depth);
1622  cz.push_back(zc);
1623  cstwd.push_back(stwd);
1624  cauxstwd.push_back(auxstwd);
1625  }
1626  }
1627 
1628  }
1629 
1630  if( subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1631 
1632  //HF
1634  ev.getByToken(tok_hf_, hfcoll);
1635 
1636  for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
1637  HcalDetId cell(j->id());
1638  const CaloCellGeometry* cellGeometry =
1639  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1640  double eta = cellGeometry->getPosition().eta () ;
1641  double phi = cellGeometry->getPosition().phi () ;
1642  double zc = cellGeometry->getPosition().z ();
1643  int sub = cell.subdet();
1644  int depth = cell.depth();
1645  int inteta = cell.ieta();
1646  if(inteta > 0) inteta -= 1;
1647  int intphi = cell.iphi()-1;
1648  double en = j->energy();
1649  double t = j->time();
1650  int stwd = j->flags();
1651  int auxstwd = j->aux();
1652 
1653  int serivityLevel = hcalSevLvl( (CaloRecHit*) &*j );
1654  if( cell.subdet()==HcalForward ){
1655  hcalHFSevLvlVec.push_back(serivityLevel);
1656  }
1657 
1658  if((iz > 0 && eta > 0.) || (iz < 0 && eta <0.) || iz == 0) {
1659 
1660  csub.push_back(sub);
1661  cen.push_back(en);
1662  ceta.push_back(eta);
1663  cphi.push_back(phi);
1664  ctime.push_back(t);
1665  cieta.push_back(inteta);
1666  ciphi.push_back(intphi);
1667  cdepth.push_back(depth);
1668  cz.push_back(zc);
1669  cstwd.push_back(stwd);
1670  cauxstwd.push_back(auxstwd);
1671  }
1672  }
1673  }
1674 
1675  //HO
1676  if( subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1677 
1679  ev.getByToken(tok_ho_, hocoll);
1680 
1681  for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
1682  HcalDetId cell(j->id());
1683  const CaloCellGeometry* cellGeometry =
1684  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1685  double eta = cellGeometry->getPosition().eta () ;
1686  double phi = cellGeometry->getPosition().phi () ;
1687  double zc = cellGeometry->getPosition().z ();
1688  int sub = cell.subdet();
1689  int depth = cell.depth();
1690  int inteta = cell.ieta();
1691  if(inteta > 0) inteta -= 1;
1692  int intphi = cell.iphi()-1;
1693  double t = j->time();
1694  double en = j->energy();
1695  int stwd = j->flags();
1696  int auxstwd = j->aux();
1697 
1698  int serivityLevel = hcalSevLvl( (CaloRecHit*) &*j );
1699  if( cell.subdet()==HcalOuter ){
1700  hcalHOSevLvlVec.push_back(serivityLevel);
1701  }
1702 
1703  if((iz > 0 && eta > 0.) || (iz < 0 && eta <0.) || iz == 0) {
1704  csub.push_back(sub);
1705  cen.push_back(en);
1706  ceta.push_back(eta);
1707  cphi.push_back(phi);
1708  ctime.push_back(t);
1709  cieta.push_back(inteta);
1710  ciphi.push_back(intphi);
1711  cdepth.push_back(depth);
1712  cz.push_back(zc);
1713  cstwd.push_back(stwd);
1714  cauxstwd.push_back(auxstwd);
1715  }
1716  }
1717  }
1718 }
1719 
1720 double HcalRecHitsValidation::dR(double eta1, double phi1, double eta2, double phi2) {
1721  double PI = 3.1415926535898;
1722  double deltaphi= phi1 - phi2;
1723  if( phi2 > phi1 ) { deltaphi= phi2 - phi1;}
1724  if(deltaphi > PI) { deltaphi = 2.*PI - deltaphi;}
1725  double deltaeta = eta2 - eta1;
1726  double tmp = sqrt(deltaeta* deltaeta + deltaphi*deltaphi);
1727  return tmp;
1728 }
1729 
1730 double HcalRecHitsValidation::phi12(double phi1, double en1, double phi2, double en2) {
1731  // weighted mean value of phi1 and phi2
1732 
1733  double tmp;
1734  double PI = 3.1415926535898;
1735  double a1 = phi1; double a2 = phi2;
1736 
1737  if( a1 > 0.5*PI && a2 < 0.) a2 += 2*PI;
1738  if( a2 > 0.5*PI && a1 < 0.) a1 += 2*PI;
1739  tmp = (a1 * en1 + a2 * en2)/(en1 + en2);
1740  if(tmp > PI) tmp -= 2.*PI;
1741 
1742  return tmp;
1743 
1744 }
1745 
1746 double HcalRecHitsValidation::dPhiWsign(double phi1, double phi2) {
1747  // clockwise phi2 w.r.t phi1 means "+" phi distance
1748  // anti-clockwise phi2 w.r.t phi1 means "-" phi distance
1749 
1750  double PI = 3.1415926535898;
1751  double a1 = phi1; double a2 = phi2;
1752  double tmp = a2 - a1;
1753  if( a1*a2 < 0.) {
1754  if(a1 > 0.5 * PI) tmp += 2.*PI;
1755  if(a2 > 0.5 * PI) tmp -= 2.*PI;
1756  }
1757  return tmp;
1758 
1759 }
1760 
1762 
1763  const DetId id = hit->detid();
1764 
1765  const uint32_t recHitFlag = hit->flags();
1766  const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1767 
1768  int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
1769 
1770  return severityLevel;
1771 
1772 }
1773 
1775 
MonitorElement * emean_vs_ieta_HB2
MonitorElement * RMS_vs_ieta_HB1
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
MonitorElement * occupancy_vs_ieta_HO
const HcalSeverityLevelComputer * theHcalSevLvlComputer
MonitorElement * meSumRecHitsEnergyConeHE
MonitorElement * occupancy_seqHB1
const HcalChannelQuality * theHcalChStatus
std::vector< PCaloHit > PCaloHitContainer
MonitorElement * occupancy_map_HB1
MonitorElement * occupancy_vs_ieta_HE3
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hh_
MonitorElement * meRecHitsEnergyHF
MonitorElement * occupancy_vs_ieta_HF1
MonitorElement * RMS_vs_ieta_HF1
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:30
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:928
#define PI
double phi12(double phi1, double en1, double phi2, double en2)
const DetId & detid() const
Definition: CaloRecHit.h:20
std::vector< int > hcalHBSevLvlVec
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
MonitorElement * meNumRecHitsConeHE
MonitorElement * meSumRecHitsEnergyConeHB
MonitorElement * meTEprofileHE_Low
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< uint32_t > cstwd
MonitorElement * RMS_vs_ieta_HE3
MonitorElement * RecHit_Aux_StatusWord_HB
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
MonitorElement * RMS_vs_ieta_HE1
MonitorElement * emean_vs_ieta_HB1
edm::EDGetTokenT< HFRecHitCollection > tok_hf_
std::vector< EcalRecHit >::const_iterator const_iterator
MonitorElement * occupancy_seqHO
MonitorElement * occupancy_seqHE3
MonitorElement * meEnConeEtaProfile
double npart
Definition: HydjetWrapper.h:44
MonitorElement * occupancy_map_HE3
double dPhiWsign(double phi1, double phi2)
MonitorElement * meRecHitSimHitProfileHF
const Item * getValues(DetId fId, bool throwOnFail=true) const
MonitorElement * meEcalHcalEnergyHB
int hcalSevLvl(const CaloRecHit *hit)
MonitorElement * meRecHitSimHitProfileHFS
MonitorElement * meEnergyHcalVsEcalHE
MonitorElement * occupancy_map_HF1
MonitorElement * meEcalHcalEnergyConeHE
MonitorElement * meRecHitsEnergyHB
MonitorElement * meSumRecHitsEnergyHF
double emap_min[82][72][4][4]
T eta() const
MonitorElement * meNumRecHitsConeHO
std::vector< double > ceta
MonitorElement * occupancy_seqHF2
MonitorElement * meNumRecHitsConeHB
MonitorElement * RMS_vs_ieta_HB2
MonitorElement * meRecHitSimHitHF
MonitorElement * meRecHitsEnergyHE
MonitorElement * occupancy_map_HE1
MonitorElement * meRecHitSimHitProfileHE
MonitorElement * meSumRecHitsEnergyHB
edm::EDGetTokenT< HORecHitCollection > tok_ho_
void Fill(long long x)
MonitorElement * meEnConeEtaProfile_E
MonitorElement * meRecHitSimHitHB
MonitorElement * meRecHitSimHitHFS
MonitorElement * RMS_vs_ieta_HO
MonitorElement * meSumRecHitsEnergyConeHFS
MonitorElement * RecHit_Aux_StatusWord_HO
std::vector< double > cphi
MonitorElement * meNumEcalRecHitsConeHB
MonitorElement * meNumRecHitsThreshHB
MonitorElement * meSumRecHitsEnergyHO
std::vector< int > hcalHESevLvlVec
MonitorElement * meEnergyHcalVsEcalHB
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
MonitorElement * RecHit_StatusWord_HF
T sqrt(T t)
Definition: SSEVec.h:48
MonitorElement * RecHit_Aux_StatusWord_HE
edm::ESHandle< CaloGeometry > geometry
MonitorElement * meSumRecHitsEnergyConeHO
MonitorElement * meTEprofileHB
uint32_t flags() const
Definition: CaloRecHit.h:21
MonitorElement * map_econe_depth1
virtual void fillRecHitsTmp(int subdet_, edm::Event const &ev)
int j
Definition: DBlmapReader.cc:9
MonitorElement * map_econe_depth4
double dR(double eta1, double phi1, double eta2, double phi2)
MonitorElement * emean_vs_ieta_HE3
MonitorElement * meNumRecHitsThreshHO
MonitorElement * meEcalHcalEnergyConeHB
MonitorElement * meNumRecHitsThreshHE
MonitorElement * bookProfile(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, const char *option="s")
Definition: DQMStore.cc:1242
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
std::vector< int > hcalHFSevLvlVec
MonitorElement * occupancy_map_HE2
MonitorElement * occupancy_vs_ieta_HB2
bool isValid() const
Definition: HandleBase.h:76
MonitorElement * meSumRecHitsEnergyConeHF
MonitorElement * RecHit_Aux_StatusWord_HF
MonitorElement * meRecHitSimHitHO
MonitorElement * meRecHitSimHitProfileHB
MonitorElement * RecHit_StatusWord_HO
MonitorElement * meEnConeEtaProfile_depth1
tuple conf
Definition: dbtoconf.py:185
MonitorElement * emean_vs_ieta_HF2
virtual void analyze(edm::Event const &ev, edm::EventSetup const &c)
MonitorElement * occupancy_map_HO
MonitorElement * occupancy_vs_ieta_HE2
MonitorElement * meTEprofileHF_Low
edm::EDGetTokenT< EBRecHitCollection > tok_EB_
MonitorElement * occupancy_map_HF2
MonitorElement * RecHit_StatusWord_HF67
MonitorElement * meTEprofileHO_High
MonitorElement * meTEprofileHE
Definition: DetId.h:18
MonitorElement * occupancy_seqHB2
MonitorElement * occupancy_map_HB2
MonitorElement * RMS_vs_ieta_HE2
MonitorElement * meNumRecHitsConeHF
MonitorElement * meTEprofileHF
MonitorElement * meEnConeEtaProfile_depth3
MonitorElement * emean_vs_ieta_HE1
MonitorElement * meRecHitSimHitProfileHO
std::vector< int > hcalHOSevLvlVec
MonitorElement * emean_vs_ieta_HE2
MonitorElement * occupancy_vs_ieta_HF2
edm::EDGetTokenT< EERecHitCollection > tok_EE_
const T & get() const
Definition: EventSetup.h:55
std::vector< double > ctime
MonitorElement * map_econe_depth2
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
MonitorElement * RecHit_StatusWord_HE
MonitorElement * meEnConeEtaProfile_EH
T const * product() const
Definition: Handle.h:81
MonitorElement * meRecHitSimHitProfileHFL
std::vector< double > cz
MonitorElement * meEnConeEtaProfile_depth2
T eta() const
Definition: PV3DBase.h:76
MonitorElement * RMS_vs_ieta_HF2
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
MonitorElement * meRecHitsEnergyHO
std::vector< uint32_t > cauxstwd
MonitorElement * emean_vs_ieta_HO
MonitorElement * occupancy_seqHF1
MonitorElement * meSumRecHitsEnergyHE
MonitorElement * meRecHitSimHitHFL
MonitorElement * meSumRecHitsEnergyConeHFL
MonitorElement * meNumEcalRecHitsConeHE
MonitorElement * meEcalHcalEnergyHE
MonitorElement * emean_vs_ieta_HF1
MonitorElement * meTEprofileHO
tuple cout
Definition: gather_cfg.py:121
edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
MonitorElement * meTEprofileHB_Low
MonitorElement * occupancy_vs_ieta_HB1
uint32_t getValue() const
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:1056
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
MonitorElement * occupancy_seqHE2
MonitorElement * meRecHitSimHitHE
HcalRecHitsValidation(edm::ParameterSet const &conf)
MonitorElement * map_econe_depth3
MonitorElement * occupancy_seqHE1
MonitorElement * meTEprofileHB_High
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:639
std::vector< double > cen
MonitorElement * meEnConeEtaProfile_depth4
MonitorElement * occupancy_vs_ieta_HE1
MonitorElement * RecHit_StatusWord_HB
Definition: DDAxes.h:10