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