test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalRecHitsValidation.cc
Go to the documentation of this file.
1 /*
2  * \file EcalRecHitsValidation.cc
3  *
4  * \author C. Rovelli
5  *
6 */
7 
14 #include <string>
18 
19 using namespace cms;
20 using namespace edm;
21 using namespace std;
22 
24 
25  // ----------------------
26  HepMCLabel = ps.getParameter<std::string>("moduleLabelMC");
27  hitsProducer_ = ps.getParameter<std::string>("hitsProducer");
28  EBrechitCollection_ = ps.getParameter<edm::InputTag>("EBrechitCollection");
29  EErechitCollection_ = ps.getParameter<edm::InputTag>("EErechitCollection");
30  ESrechitCollection_ = ps.getParameter<edm::InputTag>("ESrechitCollection");
31  EBuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EBuncalibrechitCollection");
32  EEuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EEuncalibrechitCollection");
33  // fix for consumes
34  HepMCLabel_Token_ = consumes<HepMCProduct>(ps.getParameter<std::string>("moduleLabelMC"));
35  EBrechitCollection_Token_ = consumes<EBRecHitCollection>(ps.getParameter<edm::InputTag>("EBrechitCollection"));
36  EErechitCollection_Token_ = consumes<EERecHitCollection>(ps.getParameter<edm::InputTag>("EErechitCollection"));
37  ESrechitCollection_Token_ = consumes<ESRecHitCollection>(ps.getParameter<edm::InputTag>("ESrechitCollection"));
38  EBuncalibrechitCollection_Token_ = consumes<EBUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EBuncalibrechitCollection"));
39  EEuncalibrechitCollection_Token_ = consumes<EEUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EEuncalibrechitCollection"));
40  EBHits_Token_ = consumes<CrossingFrame<PCaloHit> >(edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEB")));
41  EEHits_Token_ = consumes<CrossingFrame<PCaloHit> >(edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEE")));
42  ESHits_Token_ = consumes<CrossingFrame<PCaloHit> >(edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsES")));
43 
44  // ----------------------
45  // DQM ROOT output
46  outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
47 
48  if ( outputFile_.size() != 0 ) {
49  LogInfo("OutputInfo") << " Ecal RecHits Task histograms will be saved to '" << outputFile_.c_str() << "'";
50  } else {
51  LogInfo("OutputInfo") << " Ecal RecHits Task histograms will NOT be saved";
52  }
53 
54  // ----------------------
55  // verbosity switch
56  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
57 
58  // ----------------------
59  meGunEnergy_ = 0;
60  meGunEta_ = 0;
61  meGunPhi_ = 0;
62  meEBRecHitSimHitRatio_ = 0;
63  meEERecHitSimHitRatio_ = 0;
64  meESRecHitSimHitRatio_ = 0;
65  meEBRecHitSimHitRatio1011_ = 0;
66  meEERecHitSimHitRatio1011_ = 0;
67  meEBRecHitSimHitRatio12_ = 0;
68  meEERecHitSimHitRatio12_ = 0;
69  meEBRecHitSimHitRatio13_ = 0;
70  meEERecHitSimHitRatio13_ = 0;
71  meEBRecHitSimHitRatioGt35_ = 0;
72  meEERecHitSimHitRatioGt35_ = 0;
73  meEBUnRecHitSimHitRatio_ = 0;
74  meEEUnRecHitSimHitRatio_ = 0;
75  meEBUnRecHitSimHitRatioGt35_ = 0;
76  meEEUnRecHitSimHitRatioGt35_ = 0;
77  meEBe5x5_ = 0;
78  meEBe5x5OverSimHits_ = 0;
79  meEBe5x5OverGun_ = 0;
80  meEEe5x5_ = 0;
81  meEEe5x5OverSimHits_ = 0;
82  meEEe5x5OverGun_ = 0;
83 
84  meEBRecHitLog10Energy_ = 0;
85  meEERecHitLog10Energy_ = 0;
86  meESRecHitLog10Energy_ = 0;
87  meEBRecHitLog10EnergyContr_ = 0;
88  meEERecHitLog10EnergyContr_ = 0;
89  meESRecHitLog10EnergyContr_ = 0;
90  meEBRecHitLog10Energy5x5Contr_ = 0;
91  meEERecHitLog10Energy5x5Contr_ = 0;
92 
93  meEBRecHitsOccupancyFlag5_6_ = 0;
94  meEBRecHitsOccupancyFlag8_9_ = 0;
95  meEERecHitsOccupancyPlusFlag5_6_ = 0;
96  meEERecHitsOccupancyMinusFlag5_6_ = 0;
97  meEERecHitsOccupancyPlusFlag8_9_ = 0;
98  meEERecHitsOccupancyMinusFlag8_9_ = 0;
99 
100  meEBRecHitFlags_ = 0;
101  meEBRecHitSimHitvsSimHitFlag5_6_ = 0;
102  meEBRecHitSimHitFlag6_ = 0;
103  meEBRecHitSimHitFlag7_ = 0;
104  meEB5x5RecHitSimHitvsSimHitFlag8_ = 0;
105 
106  meEERecHitFlags_ = 0;
107  meEERecHitSimHitvsSimHitFlag5_6_ = 0;
108  meEERecHitSimHitFlag6_ = 0;
109  meEERecHitSimHitFlag7_ = 0;
110 
111 }
112 
114 
115 }
116 
118 
120 
121  ibooker.setCurrentFolder("EcalRecHitsV/EcalRecHitsTask");
122 
123  histo = "EcalRecHitsTask Gun Momentum";
124  meGunEnergy_ = ibooker.book1D(histo.c_str(), histo.c_str(), 100, 0., 1000.);
125 
126  histo = "EcalRecHitsTask Gun Eta";
127  meGunEta_ = ibooker.book1D(histo.c_str(), histo.c_str(), 700, -3.5, 3.5);
128 
129  histo = "EcalRecHitsTask Gun Phi";
130  meGunPhi_ = ibooker.book1D(histo.c_str(), histo.c_str(), 360, 0., 360.);
131 
132  histo = "EcalRecHitsTask Barrel RecSimHit Ratio";
133  meEBRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
134 
135  histo = "EcalRecHitsTask Endcap RecSimHit Ratio";
136  meEERecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
137 
138  histo = "EcalRecHitsTask Preshower RecSimHit Ratio";
139  meESRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
140 
141  histo = "EcalRecHitsTask Barrel RecSimHit Ratio gt 3p5 GeV";
142  meEBRecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
143 
144  histo = "EcalRecHitsTask Endcap RecSimHit Ratio gt 3p5 GeV";
145  meEERecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
146 
147  histo = "EcalRecHitsTask Barrel Unc RecSimHit Ratio";
148  meEBUnRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
149 
150  histo = "EcalRecHitsTask Endcap Unc RecSimHit Ratio";
151  meEEUnRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
152 
153  histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=10 11";
154  meEBRecHitSimHitRatio1011_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
155 
156  histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=10 11";
157  meEERecHitSimHitRatio1011_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
158 
159  histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=12";
160  meEBRecHitSimHitRatio12_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
161 
162  histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=12";
163  meEERecHitSimHitRatio12_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
164 
165  histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=13";
166  meEBRecHitSimHitRatio13_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
167 
168  histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=13";
169  meEERecHitSimHitRatio13_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
170 
171  histo = "EcalRecHitsTask Barrel Unc RecSimHit Ratio gt 3p5 GeV";
172  meEBUnRecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
173 
174  histo = "EcalRecHitsTask Endcap Unc RecSimHit Ratio gt 3p5 GeV";
175  meEEUnRecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
176 
177  histo = "EcalRecHitsTask Barrel Rec E5x5";
178  meEBe5x5_ = ibooker.book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
179 
180  histo = "EcalRecHitsTask Barrel Rec E5x5 over Sim E5x5";
181  meEBe5x5OverSimHits_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
182 
183  histo = "EcalRecHitsTask Barrel Rec E5x5 over gun energy";
184  meEBe5x5OverGun_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
185 
186  histo = "EcalRecHitsTask Endcap Rec E5x5";
187  meEEe5x5_ = ibooker.book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
188 
189  histo = "EcalRecHitsTask Endcap Rec E5x5 over Sim E5x5";
190  meEEe5x5OverSimHits_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
191 
192  histo = "EcalRecHitsTask Endcap Rec E5x5 over gun energy";
193  meEEe5x5OverGun_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
194 
195  meEBRecHitLog10Energy_ = ibooker.book1D( "EcalRecHitsTask Barrel Log10 Energy", "EcalRecHitsTask Barrel Log10 Energy", 90, -5., 4. );
196  meEERecHitLog10Energy_ = ibooker.book1D( "EcalRecHitsTask Endcap Log10 Energy", "EcalRecHitsTask Endcap Log10 Energy", 90, -5., 4. );
197  meESRecHitLog10Energy_ = ibooker.book1D( "EcalRecHitsTask Preshower Log10 Energy", "EcalRecHitsTask Preshower Log10 Energy", 90, -5., 4. );
198  meEBRecHitLog10EnergyContr_ = ibooker.bookProfile( "EcalRecHits Barrel Log10En vs Hit Contribution", "EcalRecHits Barrel Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
199  meEERecHitLog10EnergyContr_ = ibooker.bookProfile( "EcalRecHits Endcap Log10En vs Hit Contribution", "EcalRecHits Endcap Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
200  meESRecHitLog10EnergyContr_ = ibooker.bookProfile( "EcalRecHits Preshower Log10En vs Hit Contribution", "EcalRecHits Preshower Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
201  meEBRecHitLog10Energy5x5Contr_ = ibooker.bookProfile( "EcalRecHits Barrel Log10En5x5 vs Hit Contribution", "EcalRecHits Barrel Log10En5x5 vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
202  meEERecHitLog10Energy5x5Contr_ = ibooker.bookProfile( "EcalRecHits Endcap Log10En5x5 vs Hit Contribution", "EcalRecHits Endcap Log10En5x5 vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
203 
204 
205  histo = "EB Occupancy Flag=5 6";
206  meEBRecHitsOccupancyFlag5_6_ = ibooker.book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
207  histo = "EB Occupancy Flag=8 9";
208  meEBRecHitsOccupancyFlag8_9_ = ibooker.book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
209 
210  histo = "EE+ Occupancy Flag=5 6";
211  meEERecHitsOccupancyPlusFlag5_6_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
212  histo = "EE- Occupancy Flag=5 6";
213  meEERecHitsOccupancyMinusFlag5_6_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
214  histo = "EE+ Occupancy Flag=8 9";
215  meEERecHitsOccupancyPlusFlag8_9_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
216  histo = "EE- Occupancy Flag=8 9";
217  meEERecHitsOccupancyMinusFlag8_9_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
218 
219 
220  histo = "EcalRecHitsTask Barrel Reco Flags";
221  meEBRecHitFlags_ = ibooker.book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);
222  histo = "EcalRecHitsTask Endcap Reco Flags";
223  meEERecHitFlags_ = ibooker.book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);
224  histo = "EcalRecHitsTask Barrel RecSimHit Ratio vs SimHit Flag=5 6";
225  meEBRecHitSimHitvsSimHitFlag5_6_ = ibooker.book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );
226  histo = "EcalRecHitsTask Endcap RecSimHit Ratio vs SimHit Flag=5 6";
227  meEERecHitSimHitvsSimHitFlag5_6_ = ibooker.book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );
228  histo = "EcalRecHitsTask Barrel RecSimHit Ratio Flag=6";
229  meEBRecHitSimHitFlag6_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
230  histo = "EcalRecHitsTask Endcap RecSimHit Ratio Flag=6";
231  meEERecHitSimHitFlag6_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
232  histo = "EcalRecHitsTask Barrel RecSimHit Ratio Flag=7";
233  meEBRecHitSimHitFlag7_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
234  histo = "EcalRecHitsTask Endcap RecSimHit Ratio Flag=7";
235  meEERecHitSimHitFlag7_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
236  histo = "EcalRecHitsTask Barrel 5x5 RecSimHit Ratio vs SimHit Flag=8";
237  meEB5x5RecHitSimHitvsSimHitFlag8_ = ibooker.book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );
238 
239 }
240 
241 
243 
244  //Temporary stuff
245 
246 
247 
248  LogInfo("EcalRecHitsTask, EventInfo: ") << " Run = " << e.id().run() << " Event = " << e.id().event();
249 
250  // ADC -> GeV Scale
252  c.get<EcalADCToGeVConstantRcd>().get(pAgc);
253  const EcalADCToGeVConstant* agc = pAgc.product();
254  const double barrelADCtoGeV_ = agc->getEBValue();
255  const double endcapADCtoGeV_ = agc->getEEValue();
256 
257  Handle<HepMCProduct> MCEvt;
258  bool skipMC = false;
259  e.getByToken(HepMCLabel_Token_, MCEvt);
260  if (!MCEvt.isValid()) { skipMC = true; }
261 
262  edm::Handle<CrossingFrame<PCaloHit> > crossingFrame;
263 
264  bool skipBarrel = false;
265  const EBUncalibratedRecHitCollection *EBUncalibRecHit =0;
266  Handle< EBUncalibratedRecHitCollection > EcalUncalibRecHitEB;
267  e.getByToken( EBuncalibrechitCollection_Token_, EcalUncalibRecHitEB);
268  if (EcalUncalibRecHitEB.isValid()) {
269  EBUncalibRecHit = EcalUncalibRecHitEB.product() ;
270  } else {
271  skipBarrel = true;
272  }
273 
274  bool skipEndcap = false;
275  const EEUncalibratedRecHitCollection *EEUncalibRecHit = 0;
276  Handle< EEUncalibratedRecHitCollection > EcalUncalibRecHitEE;
277  e.getByToken( EEuncalibrechitCollection_Token_, EcalUncalibRecHitEE);
278  if (EcalUncalibRecHitEE.isValid()){
279  EEUncalibRecHit = EcalUncalibRecHitEE.product () ;
280  } else {
281  skipEndcap = true;
282  }
283 
284  const EBRecHitCollection *EBRecHit = 0;
285  Handle<EBRecHitCollection> EcalRecHitEB;
286  e.getByToken( EBrechitCollection_Token_, EcalRecHitEB);
287  if (EcalRecHitEB.isValid()){
288  EBRecHit = EcalRecHitEB.product();
289  } else {
290  skipBarrel = true;
291  }
292 
293  const EERecHitCollection *EERecHit = 0;
294  Handle<EERecHitCollection> EcalRecHitEE;
295  e.getByToken( EErechitCollection_Token_, EcalRecHitEE);
296  if (EcalRecHitEE.isValid()){
297  EERecHit = EcalRecHitEE.product ();
298  } else {
299  skipEndcap = true;
300  }
301 
302  bool skipPreshower = false;
303  const ESRecHitCollection *ESRecHit = 0;
304  Handle<ESRecHitCollection> EcalRecHitES;
305  e.getByToken( ESrechitCollection_Token_, EcalRecHitES);
306  if (EcalRecHitES.isValid()) {
307  ESRecHit = EcalRecHitES.product ();
308  } else {
309  skipPreshower = true;
310  }
311 
312 
313  // ----------------------
314  // gun
315  double eGun = 0.;
316  if ( ! skipMC ) {
317  for ( HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin(); p != MCEvt->GetEvent()->particles_end(); ++p ) {
318  double htheta = (*p)->momentum().theta();
319  double heta = -99999.;
320  if( tan(htheta * 0.5) > 0 ) {
321  heta = -log(tan(htheta * 0.5));
322  }
323  double hphi = (*p)->momentum().phi();
324  hphi = (hphi>=0) ? hphi : hphi+2*M_PI;
325  hphi = hphi / M_PI * 180.;
326 
327  LogDebug("EventInfo") << "EcalRecHitsTask: Particle gun type form MC = " << abs((*p)->pdg_id())
328  << "\n" << "Energy = "<< (*p)->momentum().e()
329  << "\n" << "Eta = " << heta
330  << "\n" << "Phi = " << hphi;
331 
332  if ( (*p)->momentum().e() > eGun ) eGun = (*p)->momentum().e();
333 
334  if (meGunEnergy_) meGunEnergy_->Fill((*p)->momentum().e());
335  if (meGunEta_) meGunEta_ ->Fill(heta);
336  if (meGunPhi_) meGunPhi_ ->Fill(hphi);
337  }
338  }
339 
340  // -------------------------------------------------------------------
341  // BARREL
342 
343  if ( ! skipBarrel) {
344 
345  // 1) loop over simHits
346  e.getByToken(EBHits_Token_,crossingFrame);
347  std::auto_ptr<MixCollection<PCaloHit> >
348  barrelHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
349 
350  MapType ebSimMap;
351  MapType ebRecMap;
352  const int ebcSize = 90;
353  double ebcontr[ebcSize];
354  double ebcontr25[ebcSize];
355  for( int i=0; i<ebcSize; i++ ) { ebcontr[i] = 0.0; ebcontr25[i] = 0.0; }
356  double ebtotal = 0.;
357 
358  for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin (); hitItr != barrelHits->end (); ++hitItr) {
359  EBDetId ebid = EBDetId(hitItr->id());
360 
361  LogDebug("SimHitInfo, barrel")
362  << "CaloHit " << hitItr->getName() << " DetID = " << hitItr->id() << "\n"
363  << "Energy = " << hitItr->energy() << " Time = " << hitItr->time() << "\n"
364  << "EBDetId = " << ebid.ieta() << " " << ebid.iphi();
365 
366  uint32_t crystid = ebid.rawId();
367  ebSimMap[crystid] += hitItr->energy();
368  }
369 
370 
371 
372  // 2) loop over RecHits
373  for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin(); uncalibRecHit != EBUncalibRecHit->end() ; ++uncalibRecHit) {
374  EBDetId EBid = EBDetId(uncalibRecHit->id());
375 
376  // Find corresponding recHit
377  EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
378  if( myRecHit == EBRecHit->end() ) continue;
379  ebRecMap[EBid.rawId()] += myRecHit->energy();
380 
381  // Fill log10(Energy) stuff...
382  ebtotal += myRecHit->energy();
383  if( myRecHit->energy() > 0 ) {
384  if( meEBRecHitLog10Energy_ ) meEBRecHitLog10Energy_->Fill( log10( myRecHit->energy() ) );
385  int log10i = int( ( log10( myRecHit->energy() ) + 5. ) * 10. );
386  if( log10i >=0 and log10i < ebcSize ) ebcontr[ log10i ] += myRecHit->energy();
387  }
388 
389  // comparison Rec/Sim hit
390  if ( ebSimMap[EBid.rawId()] != 0. ) {
391  double uncEnergy = uncalibRecHit->amplitude()*barrelADCtoGeV_;
392  if (meEBUnRecHitSimHitRatio_) {meEBUnRecHitSimHitRatio_ ->Fill(uncEnergy/ebSimMap[EBid.rawId()]);}
393  if (meEBUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBUnRecHitSimHitRatioGt35_->Fill(uncEnergy/ebSimMap[EBid.rawId()]);}
394  }
395 
396  if (myRecHit != EBRecHit->end()) {
397  if ( ebSimMap[EBid.rawId()] != 0. ) {
398  if (meEBRecHitSimHitRatio_) {meEBRecHitSimHitRatio_ ->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]);}
399  if (meEBRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBRecHitSimHitRatioGt35_->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]);}
400  uint16_t sc = 0;
402  c.get<EcalChannelStatusRcd>().get(pEcs);
403  const EcalChannelStatus* ecs = 0;
404  if( pEcs.isValid() ) ecs = pEcs.product();
405  if( ecs != 0 ) {
406  EcalChannelStatusMap::const_iterator csmi = ecs->find(EBid.rawId());
407  EcalChannelStatusCode csc = 0;
408  if( csmi != ecs->end() ) csc = *csmi;
409  sc = csc.getStatusCode();
410  }
411 
412  if( meEBRecHitSimHitRatio1011_ != 0 &&
413  ( sc == 10 || sc == 11 ) ) { meEBRecHitSimHitRatio1011_->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]); }
414  if( meEBRecHitSimHitRatio12_ != 0 && sc == 12 ) { meEBRecHitSimHitRatio12_->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]); }
415 
417  c.get<IdealGeometryRecord>().get(pttMap);
418  const EcalTrigTowerConstituentsMap* ttMap = 0;
419  if( pttMap.isValid() ) ttMap = pttMap.product();
420  double ttSimEnergy = 0;
421  if( ttMap != 0 ) {
422  EcalTrigTowerDetId ttDetId = EBid.tower();
423  std::vector<DetId> vid = ttMap->constituentsOf( ttDetId );
424  for( std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); dit++ ) {
425  EBDetId ttEBid = EBDetId(*dit);
426  ttSimEnergy += ebSimMap[ttEBid.rawId()];
427  }
428  if( vid.size() != 0 ) ttSimEnergy = ttSimEnergy / vid.size();
429  }
430  if( meEBRecHitSimHitRatio13_ != 0 && sc == 13 && ttSimEnergy != 0 )
431  meEBRecHitSimHitRatio13_->Fill(myRecHit->energy()/ttSimEnergy);
432 
433  int flag = myRecHit->recoFlag();
434  if( meEBRecHitFlags_ != 0 ) meEBRecHitFlags_->Fill( flag );
435  if( meEBRecHitSimHitvsSimHitFlag5_6_ && ( flag == EcalRecHit::kSaturated || flag == EcalRecHit::kLeadingEdgeRecovered ))
436  meEBRecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/ebSimMap[EBid.rawId()], ebSimMap[EBid.rawId()] );
437  if( meEBRecHitSimHitFlag6_ && ( flag == EcalRecHit::kLeadingEdgeRecovered ))
438  meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.rawId()] );
439  if( meEBRecHitSimHitFlag7_ && ( flag == EcalRecHit::kNeighboursRecovered ))
440  meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.rawId()] );
441  if( meEB5x5RecHitSimHitvsSimHitFlag8_ && ( flag == EcalRecHit::kTowerRecovered ) && ttSimEnergy != 0 )
442  meEB5x5RecHitSimHitvsSimHitFlag8_->Fill( myRecHit->energy()/ttSimEnergy, ttSimEnergy );
443 
444  if (meEBRecHitsOccupancyFlag5_6_ && ( (flag==EcalRecHit::kSaturated) || (flag==EcalRecHit::kLeadingEdgeRecovered) ) )
445  meEBRecHitsOccupancyFlag5_6_ -> Fill(EBid.ieta(), EBid.iphi());
446  if (meEBRecHitsOccupancyFlag8_9_ && ( (flag==EcalRecHit::kTowerRecovered) || (flag==EcalRecHit::kDead) ) )
447  meEBRecHitsOccupancyFlag8_9_ -> Fill(EBid.ieta(), EBid.iphi());
448 
449  }
450  }
451  else
452  continue;
453  } // loop over the UncalibratedRecHitCollection
454 
455  // RecHits matrix
456  uint32_t ebcenterid = getUnitWithMaxEnergy(ebRecMap);
457  EBDetId myEBid(ebcenterid);
458  int bx = myEBid.ietaAbs();
459  int by = myEBid.iphi();
460  int bz = myEBid.zside();
461  findBarrelMatrix(5,5,bx,by,bz,ebRecMap);
462  double e5x5rec = 0.;
463  double e5x5sim = 0.;
464  for ( unsigned int i = 0; i < crystalMatrix.size(); i++ ) {
465  e5x5rec += ebRecMap[crystalMatrix[i]];
466  e5x5sim += ebSimMap[crystalMatrix[i]];
467  if( ebRecMap[crystalMatrix[i]] > 0 ) {
468  int log10i25 = int( ( log10( ebRecMap[crystalMatrix[i]] ) + 5. ) * 10. );
469  if( log10i25 >=0 && log10i25 < ebcSize ) ebcontr25[ log10i25 ] += ebRecMap[crystalMatrix[i]];
470  }
471  }
472 
473  if( meEBe5x5_ ) meEBe5x5_->Fill(e5x5rec);
474  if ( e5x5sim > 0. && meEBe5x5OverSimHits_ ) meEBe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
475  if ( eGun > 0. && meEBe5x5OverGun_ ) meEBe5x5OverGun_->Fill(e5x5rec/eGun);
476 
477  if( meEBRecHitLog10EnergyContr_ && ebtotal != 0 ) {
478  for( int i=0; i<ebcSize; i++ ) {
479  meEBRecHitLog10EnergyContr_->Fill( -5.+(float(i)+0.5)/10., ebcontr[i]/ebtotal );
480  }
481  }
482 
483  if( meEBRecHitLog10Energy5x5Contr_ && e5x5rec != 0 ) {
484  for( int i=0; i<ebcSize; i++ ) {
485  meEBRecHitLog10Energy5x5Contr_->Fill( -5.+(float(i)+0.5)/10., ebcontr25[i]/e5x5rec );
486  }
487  }
488  }
489 
490  // -------------------------------------------------------------------
491  // ENDCAP
492 
493  if ( ! skipEndcap ) {
494 
495  // 1) loop over simHits
496  e.getByToken(EEHits_Token_,crossingFrame);
497  std::auto_ptr<MixCollection<PCaloHit> >
498  endcapHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
499 
500  MapType eeSimMap;
501  MapType eeRecMap;
502  const int eecSize = 90;
503  double eecontr[eecSize];
504  double eecontr25[eecSize];
505  for( int i=0; i<eecSize; i++ ) { eecontr[i] = 0.0; eecontr25[i] = 0.0; }
506  double eetotal = 0.;
507 
508  for (MixCollection<PCaloHit>::MixItr hitItr = endcapHits->begin(); hitItr != endcapHits->end(); ++hitItr) {
509  EEDetId eeid = EEDetId(hitItr->id()) ;
510 
511  LogDebug("Endcap, HitInfo")
512  <<" CaloHit " << hitItr->getName() << " DetID = " << hitItr->id() << "\n"
513  << "Energy = " << hitItr->energy() << " Time = " << hitItr->time() << "\n"
514  << "EEDetId side " << eeid.zside() << " = " << eeid.ix() << " " << eeid.iy();
515 
516  uint32_t crystid = eeid.rawId();
517  eeSimMap[crystid] += hitItr->energy();
518  }
519 
520 
521 
522  // 2) loop over RecHits
523  for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EEUncalibRecHit->begin(); uncalibRecHit != EEUncalibRecHit->end(); ++uncalibRecHit) {
524  EEDetId EEid = EEDetId(uncalibRecHit->id());
525 
526  // Find corresponding recHit
527  EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
528  if( myRecHit == EERecHit->end() ) continue;
529  eeRecMap[EEid.rawId()] += myRecHit->energy();
530 
531  // Fill log10(Energy) stuff...
532  eetotal += myRecHit->energy();
533  if( myRecHit->energy() > 0 ) {
534  if( meEERecHitLog10Energy_ ) meEERecHitLog10Energy_->Fill( log10( myRecHit->energy() ) );
535  int log10i = int( ( log10( myRecHit->energy() ) + 5. ) * 10. );
536  if( log10i >=0 and log10i < eecSize ) eecontr[ log10i ] += myRecHit->energy();
537  }
538 
539  // comparison Rec/Sim hit
540  if ( eeSimMap[EEid.rawId()] != 0. ) {
541  double uncEnergy = uncalibRecHit->amplitude()*endcapADCtoGeV_;
542  if (meEEUnRecHitSimHitRatio_) {meEEUnRecHitSimHitRatio_ ->Fill(uncEnergy/eeSimMap[EEid.rawId()]);}
543  if (meEEUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEEUnRecHitSimHitRatioGt35_->Fill(uncEnergy/eeSimMap[EEid.rawId()]);}
544  }
545 
546  if (myRecHit != EERecHit->end()) {
547  if ( eeSimMap[EEid.rawId()] != 0. ) {
548  if (meEERecHitSimHitRatio_) {meEERecHitSimHitRatio_ ->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
549  if (meEERecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEERecHitSimHitRatioGt35_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
550 
552  c.get<EcalChannelStatusRcd>().get(pEcs);
553  const EcalChannelStatus* ecs = 0;
554  if( pEcs.isValid() ) ecs = pEcs.product();
555  if( ecs != 0 ) {
556  EcalChannelStatusMap::const_iterator csmi = ecs->find(EEid.rawId());
557  EcalChannelStatusCode csc = 0;
558  if( csmi != ecs->end() ) csc = *csmi;
559  uint16_t sc = csc.getStatusCode();
560  if( meEERecHitSimHitRatio1011_ != 0 &&
561  ( sc == 10 || sc == 11 ) ) { meEERecHitSimHitRatio1011_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
562  if( meEERecHitSimHitRatio12_ != 0 && sc == 12 ) { meEERecHitSimHitRatio12_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
563  if( meEERecHitSimHitRatio13_ != 0 && sc == 13 ) { meEERecHitSimHitRatio13_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
564  }
565 
566  int flag = myRecHit->recoFlag();
567  if( meEERecHitFlags_ != 0 ) meEERecHitFlags_->Fill( flag );
568  if( meEERecHitSimHitvsSimHitFlag5_6_ && ( flag == EcalRecHit::kSaturated || flag == EcalRecHit::kLeadingEdgeRecovered ))
569  meEERecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/eeSimMap[EEid.rawId()], eeSimMap[EEid.rawId()] );
570  if( meEERecHitSimHitFlag6_ && ( flag == EcalRecHit::kLeadingEdgeRecovered ))
571  meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.rawId()] );
572  if( meEERecHitSimHitFlag7_ && ( flag == EcalRecHit::kNeighboursRecovered ))
573  meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.rawId()] );
574 
575  if (EEid.zside() > 0) {
576  if (meEERecHitsOccupancyPlusFlag5_6_ && (( flag == EcalRecHit::kSaturated ) || ( flag == EcalRecHit::kLeadingEdgeRecovered ) ))
577  meEERecHitsOccupancyPlusFlag5_6_ ->Fill(EEid.ix(), EEid.iy());
578  if (meEERecHitsOccupancyPlusFlag8_9_ && (( flag == EcalRecHit::kTowerRecovered ) || ( flag == EcalRecHit::kDead ) ))
579  meEERecHitsOccupancyPlusFlag8_9_ ->Fill(EEid.ix(), EEid.iy());
580  }
581  if (EEid.zside() < 0) {
582  if (meEERecHitsOccupancyMinusFlag5_6_ && (( flag == EcalRecHit::kSaturated ) || ( flag == EcalRecHit::kLeadingEdgeRecovered ) ))
583  meEERecHitsOccupancyMinusFlag5_6_ ->Fill(EEid.ix(), EEid.iy());
584  if (meEERecHitsOccupancyMinusFlag8_9_ && (( flag == EcalRecHit::kTowerRecovered ) || ( flag == EcalRecHit::kDead ) ))
585  meEERecHitsOccupancyMinusFlag8_9_ ->Fill(EEid.ix(), EEid.iy());
586  }
587  }
588  }
589  else
590  continue;
591  } // loop over the UncalibratedechitCollection
592 
593  // RecHits matrix
594  uint32_t eecenterid = getUnitWithMaxEnergy(eeRecMap);
595  EEDetId myEEid(eecenterid);
596  int bx = myEEid.ix();
597  int by = myEEid.iy();
598  int bz = myEEid.zside();
599  findEndcapMatrix(5,5,bx,by,bz,eeRecMap);
600  double e5x5rec = 0.;
601  double e5x5sim = 0.;
602  for ( unsigned int i = 0; i < crystalMatrix.size(); i++ ) {
603  e5x5rec += eeRecMap[crystalMatrix[i]];
604  e5x5sim += eeSimMap[crystalMatrix[i]];
605  if( eeRecMap[crystalMatrix[i]] > 0 ) {
606  int log10i25 = int( ( log10( eeRecMap[crystalMatrix[i]] ) + 5. ) * 10. );
607  if( log10i25 >=0 && log10i25 < eecSize ) eecontr25[ log10i25 ] += eeRecMap[crystalMatrix[i]];
608  }
609  }
610 
611  if( meEEe5x5_ ) meEEe5x5_->Fill(e5x5rec);
612  if ( e5x5sim > 0. && meEEe5x5OverSimHits_ ) meEEe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
613  if ( eGun > 0. && meEEe5x5OverGun_ ) meEEe5x5OverGun_->Fill(e5x5rec/eGun);
614 
615 
616  if( meEERecHitLog10EnergyContr_ && eetotal != 0 ) {
617  for( int i=0; i<eecSize; i++ ) {
618  meEERecHitLog10EnergyContr_->Fill( -5.+(float(i)+0.5)/10., eecontr[i]/eetotal );
619  }
620  }
621 
622  if( meEERecHitLog10Energy5x5Contr_ && e5x5rec != 0 ) {
623  for( int i=0; i<eecSize; i++ ) {
624  meEERecHitLog10Energy5x5Contr_->Fill( -5.+(float(i)+0.5)/10., eecontr25[i]/e5x5rec );
625  }
626  }
627  }
628 
629  // -------------------------------------------------------------------
630  // PRESHOWER
631 
632  if ( ! skipPreshower ) {
633 
634  // 1) loop over simHits
635  e.getByToken(ESHits_Token_,crossingFrame);
636  std::auto_ptr<MixCollection<PCaloHit> >
637  preshowerHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
638 
639  MapType esSimMap;
640  const int escSize = 90;
641  double escontr[escSize];
642  for( int i=0; i<escSize; i++ ) { escontr[i] = 0.0; }
643  double estotal = 0.;
644 
645 
646  for (MixCollection<PCaloHit>::MixItr hitItr = preshowerHits->begin(); hitItr != preshowerHits->end(); ++hitItr) {
647  ESDetId esid = ESDetId(hitItr->id()) ;
648 
649  LogDebug("Preshower, HitInfo")
650  <<" CaloHit " << hitItr->getName() << " DetID = " << hitItr->id() << "\n"
651  << "Energy = " << hitItr->energy() << " Time = " << hitItr->time() << "\n"
652  << "ESDetId strip " << esid.strip() << " = " << esid.six() << " " << esid.siy();
653 
654  uint32_t crystid = esid.rawId();
655  esSimMap[crystid] += hitItr->energy();
656  }
657 
658 
659  // 2) loop over RecHits
660  for (EcalRecHitCollection::const_iterator recHit = ESRecHit->begin(); recHit != ESRecHit->end(); ++recHit) {
661  ESDetId ESid = ESDetId(recHit->id());
662  if ( esSimMap[ESid.rawId()] != 0. ) {
663 
664  // Fill log10(Energy) stuff...
665  estotal += recHit->energy();
666  if( recHit->energy() > 0 ) {
667  if( meESRecHitLog10Energy_ ) meESRecHitLog10Energy_->Fill( log10( recHit->energy() ) );
668  int log10i = int( ( log10( recHit->energy() ) + 5. ) * 10. );
669  if( log10i >=0 and log10i < escSize ) escontr[ log10i ] += recHit->energy();
670  }
671 
672  if (meESRecHitSimHitRatio_) {
673  meESRecHitSimHitRatio_ ->Fill(recHit->energy()/esSimMap[ESid.rawId()]);
674  }
675  }
676  else
677  continue;
678  } // loop over the RechitCollection
679 
680  if( meESRecHitLog10EnergyContr_ && estotal != 0 ) {
681  for( int i=0; i<escSize; i++ ) {
682  meESRecHitLog10EnergyContr_->Fill( -5.+(float(i)+0.5)/10., escontr[i]/estotal );
683  }
684  }
685 
686 
687  }
688 
689 }
690 
691 
693 
694  //look for max
695  uint32_t unitWithMaxEnergy = 0;
696  float maxEnergy = 0.;
697 
698  MapType::iterator iter;
699  for (iter = themap.begin(); iter != themap.end(); iter++) {
700 
701  if (maxEnergy < (*iter).second) {
702  maxEnergy = (*iter).second;
703  unitWithMaxEnergy = (*iter).first;
704  }
705  }
706 
707  return unitWithMaxEnergy;
708 }
709 
710 void EcalRecHitsValidation::findBarrelMatrix(int nCellInEta, int nCellInPhi,
711  int CentralEta, int CentralPhi,int CentralZ,
712  MapType& themap) {
713 
714  int goBackInEta = nCellInEta/2;
715  int goBackInPhi = nCellInPhi/2;
716  int matrixSize = nCellInEta*nCellInPhi;
717  crystalMatrix.clear();
718  crystalMatrix.resize(matrixSize);
719 
720  int startEta = CentralZ*CentralEta - goBackInEta;
721  int startPhi = CentralPhi - goBackInPhi;
722 
723  int i = 0 ;
724  for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
725  for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
726  uint32_t index;
727  if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
728  if (iphi< 1) { index = EBDetId(ieta,iphi+360).rawId(); }
729  else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
730  else { index = EBDetId(ieta,iphi).rawId(); }
731  crystalMatrix[i++] = index;
732  }
733  }
734 
735 }
736 
737 void EcalRecHitsValidation::findEndcapMatrix(int nCellInX, int nCellInY,
738  int CentralX, int CentralY,int CentralZ,
739  MapType& themap) {
740  int goBackInX = nCellInX/2;
741  int goBackInY = nCellInY/2;
742  crystalMatrix.clear();
743 
744  int startX = CentralX - goBackInX;
745  int startY = CentralY - goBackInY;
746 
747  for ( int ix = startX; ix < startX+nCellInX; ix ++ ) {
748 
749  for( int iy = startY; iy < startY + nCellInY; iy++ ) {
750 
751  uint32_t index ;
752 
753  if(EEDetId::validDetId(ix,iy,CentralZ)) {
754  index = EEDetId(ix,iy,CentralZ).rawId();
755  }
756  else { continue; }
757  crystalMatrix.push_back(index);
758  }
759  }
760 }
#define LogDebug(id)
RunNumber_t run() const
Definition: EventID.h:39
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
int strip() const
Definition: ESDetId.h:52
int ix() const
Definition: EEDetId.h:76
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
EcalRecHitsValidation(const edm::ParameterSet &ps)
Constructor.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
void analyze(const edm::Event &e, const edm::EventSetup &c)
Analyze.
Code getStatusCode() const
return decoded status
std::vector< EcalUncalibratedRecHit >::const_iterator const_iterator
int six() const
Definition: ESDetId.h:48
uint32_t getUnitWithMaxEnergy(MapType &themap)
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
int siy() const
Definition: ESDetId.h:50
void findEndcapMatrix(int nCellInX, int nCellInY, int CentralX, int CentralY, int CentralZ, MapType &themap)
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
int zside() const
Definition: EEDetId.h:70
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int iy() const
Definition: EEDetId.h:82
EcalTrigTowerDetId tower() const
get the HCAL/trigger iphi of this crystal
Definition: EBDetId.h:59
int ieta() const
get the crystal ieta
Definition: EBDetId.h:51
bool isValid() const
Definition: HandleBase.h:75
std::vector< DetId > constituentsOf(const EcalTrigTowerDetId &id) const
Get the constituent detids for this tower id.
void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
#define M_PI
const_iterator end() const
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:276
T const * product() const
Definition: Handle.h:81
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:248
const T & get() const
Definition: EventSetup.h:56
std::vector< Item >::const_iterator const_iterator
T const * product() const
Definition: ESHandle.h:86
edm::EventID id() const
Definition: EventBase.h:59
iterator find(key_type k)
void findBarrelMatrix(int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &themap)
const_iterator find(uint32_t rawId) const
const_iterator end() const
std::map< uint32_t, float, std::less< uint32_t > > MapType
bool isValid() const
Definition: ESHandle.h:47
int ietaAbs() const
get the absolute value of the crystal ieta
Definition: EBDetId.h:49
const_iterator begin() const
Definition: Run.h:43
int zside() const
get the z-side of the crystal (1/-1)
Definition: EBDetId.h:47