27 HepMCLabel = ps.
getParameter<std::string>(
"moduleLabelMC");
28 hitsProducer_ = ps.
getParameter<std::string>(
"hitsProducer");
39 if ( outputFile_.size() != 0 ) {
40 LogInfo(
"OutputInfo") <<
" Ecal RecHits Task histograms will be saved to '" << outputFile_.c_str() <<
"'";
42 LogInfo(
"OutputInfo") <<
" Ecal RecHits Task histograms will NOT be saved";
61 if ( verbose_ ) dbe_->showDirStructure();
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;
85 meEBe5x5OverSimHits_ = 0;
88 meEEe5x5OverSimHits_ = 0;
91 meEBRecHitLog10Energy_ = 0;
92 meEERecHitLog10Energy_ = 0;
93 meESRecHitLog10Energy_ = 0;
94 meEBRecHitLog10EnergyContr_ = 0;
95 meEERecHitLog10EnergyContr_ = 0;
96 meESRecHitLog10EnergyContr_ = 0;
97 meEBRecHitLog10Energy5x5Contr_ = 0;
98 meEERecHitLog10Energy5x5Contr_ = 0;
100 meEBRecHitsOccupancyFlag5_6_ = 0;
101 meEBRecHitsOccupancyFlag8_9_ = 0;
102 meEERecHitsOccupancyPlusFlag5_6_ = 0;
103 meEERecHitsOccupancyMinusFlag5_6_ = 0;
104 meEERecHitsOccupancyPlusFlag8_9_ = 0;
105 meEERecHitsOccupancyMinusFlag8_9_ = 0;
107 meEBRecHitFlags_ = 0;
108 meEBRecHitSimHitvsSimHitFlag5_6_ = 0;
109 meEBRecHitSimHitFlag6_ = 0;
110 meEBRecHitSimHitFlag7_ = 0;
111 meEB5x5RecHitSimHitvsSimHitFlag8_ = 0;
113 meEERecHitFlags_ = 0;
114 meEERecHitSimHitvsSimHitFlag5_6_ = 0;
115 meEERecHitSimHitFlag6_ = 0;
116 meEERecHitSimHitFlag7_ = 0;
123 dbe_->setCurrentFolder(
"EcalRecHitsV/EcalRecHitsTask");
125 histo =
"EcalRecHitsTask Gun Momentum";
126 meGunEnergy_ = dbe_->book1D(histo.c_str(), histo.c_str(), 100, 0., 1000.);
128 histo =
"EcalRecHitsTask Gun Eta";
129 meGunEta_ = dbe_->book1D(histo.c_str(), histo.c_str(), 700, -3.5, 3.5);
131 histo =
"EcalRecHitsTask Gun Phi";
132 meGunPhi_ = dbe_->book1D(histo.c_str(), histo.c_str(), 360, 0., 360.);
134 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio";
135 meEBRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
137 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio";
138 meEERecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
140 histo =
"EcalRecHitsTask Preshower RecSimHit Ratio";
141 meESRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
143 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio gt 3p5 GeV";
144 meEBRecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
146 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio gt 3p5 GeV";
147 meEERecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
149 histo =
"EcalRecHitsTask Barrel Unc RecSimHit Ratio";
150 meEBUnRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
152 histo =
"EcalRecHitsTask Endcap Unc RecSimHit Ratio";
153 meEEUnRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
155 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=10 11";
156 meEBRecHitSimHitRatio1011_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
158 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=10 11";
159 meEERecHitSimHitRatio1011_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
161 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=12";
162 meEBRecHitSimHitRatio12_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
164 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=12";
165 meEERecHitSimHitRatio12_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
167 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=13";
168 meEBRecHitSimHitRatio13_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
170 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=13";
171 meEERecHitSimHitRatio13_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
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);
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);
179 histo =
"EcalRecHitsTask Barrel Rec E5x5";
180 meEBe5x5_ = dbe_->book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
182 histo =
"EcalRecHitsTask Barrel Rec E5x5 over Sim E5x5";
183 meEBe5x5OverSimHits_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
185 histo =
"EcalRecHitsTask Barrel Rec E5x5 over gun energy";
186 meEBe5x5OverGun_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
188 histo =
"EcalRecHitsTask Endcap Rec E5x5";
189 meEEe5x5_ = dbe_->book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
191 histo =
"EcalRecHitsTask Endcap Rec E5x5 over Sim E5x5";
192 meEEe5x5OverSimHits_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
194 histo =
"EcalRecHitsTask Endcap Rec E5x5 over gun energy";
195 meEEe5x5OverGun_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
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. );
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.);
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.);
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. );
246 if ( outputFile_.size() != 0 &&
dbe_ )
dbe_->
save(outputFile_);
263 LogInfo(
"EcalRecHitsTask, EventInfo: ") <<
" Run = " << e.
id().
run() <<
" Event = " << e.
id().
event();
269 const double barrelADCtoGeV_ = agc->
getEBValue();
270 const double endcapADCtoGeV_ = agc->
getEEValue();
275 if (!MCEvt.
isValid()) { skipMC =
true; }
279 bool skipBarrel =
false;
282 e.
getByLabel( EBuncalibrechitCollection_, EcalUncalibRecHitEB);
283 if (EcalUncalibRecHitEB.
isValid()) {
284 EBUncalibRecHit = EcalUncalibRecHitEB.
product() ;
289 bool skipEndcap =
false;
292 e.
getByLabel( EEuncalibrechitCollection_, EcalUncalibRecHitEE);
293 if (EcalUncalibRecHitEE.
isValid()){
294 EEUncalibRecHit = EcalUncalibRecHitEE.
product () ;
301 e.
getByLabel( EBrechitCollection_, EcalRecHitEB);
303 EBRecHit = EcalRecHitEB.
product();
310 e.
getByLabel( EErechitCollection_, EcalRecHitEE);
312 EERecHit = EcalRecHitEE.
product ();
317 bool skipPreshower =
false;
320 e.
getByLabel( ESrechitCollection_, EcalRecHitES);
322 ESRecHit = EcalRecHitES.
product ();
324 skipPreshower =
true;
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));
338 double hphi = (*p)->momentum().phi();
339 hphi = (hphi>=0) ? hphi : hphi+2*
M_PI;
340 hphi = hphi /
M_PI * 180.;
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;
347 if ( (*p)->momentum().e() > eGun ) eGun = (*p)->momentum().e();
349 if (meGunEnergy_) meGunEnergy_->Fill((*p)->momentum().e());
350 if (meGunEta_) meGunEta_ ->Fill(heta);
351 if (meGunPhi_) meGunPhi_ ->Fill(hphi);
361 const std::string barrelHitsName(hitsProducer_+
"EcalHitsEB");
362 e.
getByLabel(
"mix",barrelHitsName,crossingFrame);
363 std::auto_ptr<MixCollection<PCaloHit> >
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; }
378 <<
"CaloHit " << hitItr->getName() <<
" DetID = " << hitItr->id() <<
"\n"
379 <<
"Energy = " << hitItr->energy() <<
" Time = " << hitItr->time() <<
"\n"
380 <<
"EBDetId = " << ebid.
ieta() <<
" " << ebid.
iphi();
382 uint32_t crystid = ebid.
rawId();
383 ebSimMap[crystid] += hitItr->energy();
394 if( myRecHit == EBRecHit->
end() )
continue;
395 ebRecMap[EBid.
rawId()] += myRecHit->energy();
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();
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()]);}
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()]);}
424 if( csmi != ecs->
end() ) csc = *csmi;
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()]); }
436 double ttSimEnergy = 0;
440 for( std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); dit++ ) {
442 ttSimEnergy += ebSimMap[ttEBid.
rawId()];
444 if( vid.size() != 0 ) ttSimEnergy = ttSimEnergy / vid.size();
446 if( meEBRecHitSimHitRatio13_ != 0 && sc == 13 && ttSimEnergy != 0 )
447 meEBRecHitSimHitRatio13_->Fill(myRecHit->energy()/ttSimEnergy);
449 int flag = myRecHit->recoFlag();
450 if( meEBRecHitFlags_ != 0 ) meEBRecHitFlags_->Fill( flag );
452 meEBRecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/ebSimMap[EBid.
rawId()], ebSimMap[EBid.
rawId()] );
454 meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.
rawId()] );
456 meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.
rawId()] );
458 meEB5x5RecHitSimHitvsSimHitFlag8_->Fill( myRecHit->energy()/ttSimEnergy, ttSimEnergy );
461 meEBRecHitsOccupancyFlag5_6_ -> Fill(EBid.
ieta(), EBid.
iphi());
463 meEBRecHitsOccupancyFlag8_9_ -> Fill(EBid.
ieta(), EBid.
iphi());
472 uint32_t ebcenterid = getUnitWithMaxEnergy(ebRecMap);
475 int by = myEBid.
iphi();
476 int bz = myEBid.
zside();
477 findBarrelMatrix(5,5,bx,by,bz,ebRecMap);
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]];
489 if( meEBe5x5_ ) meEBe5x5_->Fill(e5x5rec);
490 if ( e5x5sim > 0. && meEBe5x5OverSimHits_ ) meEBe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
491 if ( eGun > 0. && meEBe5x5OverGun_ ) meEBe5x5OverGun_->Fill(e5x5rec/eGun);
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 );
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 );
509 if ( ! skipEndcap ) {
512 const std::string endcapHitsName(hitsProducer_+
"EcalHitsEE");
513 e.
getByLabel(
"mix",endcapHitsName,crossingFrame);
514 std::auto_ptr<MixCollection<PCaloHit> >
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; }
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();
533 uint32_t crystid = eeid.
rawId();
534 eeSimMap[crystid] += hitItr->energy();
545 if( myRecHit == EERecHit->
end() )
continue;
546 eeRecMap[EEid.
rawId()] += myRecHit->energy();
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();
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()]);}
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()]); }
575 if( csmi != ecs->
end() ) csc = *csmi;
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()]); }
583 int flag = myRecHit->recoFlag();
584 if( meEERecHitFlags_ != 0 ) meEERecHitFlags_->Fill( flag );
586 meEERecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/eeSimMap[EEid.
rawId()], eeSimMap[EEid.
rawId()] );
588 meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.
rawId()] );
590 meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.
rawId()] );
592 if (EEid.
zside() > 0) {
594 meEERecHitsOccupancyPlusFlag5_6_ ->Fill(EEid.
ix(), EEid.
iy());
596 meEERecHitsOccupancyPlusFlag8_9_ ->Fill(EEid.
ix(), EEid.
iy());
598 if (EEid.
zside() < 0) {
600 meEERecHitsOccupancyMinusFlag5_6_ ->Fill(EEid.
ix(), EEid.
iy());
602 meEERecHitsOccupancyMinusFlag8_9_ ->Fill(EEid.
ix(), EEid.
iy());
611 uint32_t eecenterid = getUnitWithMaxEnergy(eeRecMap);
613 int bx = myEEid.
ix();
614 int by = myEEid.
iy();
615 int bz = myEEid.
zside();
616 findEndcapMatrix(5,5,bx,by,bz,eeRecMap);
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]];
628 if( meEEe5x5_ ) meEEe5x5_->Fill(e5x5rec);
629 if ( e5x5sim > 0. && meEEe5x5OverSimHits_ ) meEEe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
630 if ( eGun > 0. && meEEe5x5OverGun_ ) meEEe5x5OverGun_->Fill(e5x5rec/eGun);
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 );
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 );
649 if ( ! skipPreshower ) {
652 const std::string preshowerHitsName(hitsProducer_+
"EcalHitsES");
653 e.
getByLabel(
"mix",preshowerHitsName,crossingFrame);
654 std::auto_ptr<MixCollection<PCaloHit> >
658 const int escSize = 90;
659 double escontr[escSize];
660 for(
int i=0;
i<escSize;
i++ ) { escontr[
i] = 0.0; }
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();
672 uint32_t crystid = esid.
rawId();
673 esSimMap[crystid] += hitItr->energy();
680 if ( esSimMap[ESid.
rawId()] != 0. ) {
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();
690 if (meESRecHitSimHitRatio_) {
691 meESRecHitSimHitRatio_ ->Fill(recHit->energy()/esSimMap[ESid.
rawId()]);
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 );
713 uint32_t unitWithMaxEnergy = 0;
714 float maxEnergy = 0.;
716 MapType::iterator iter;
717 for (iter = themap.begin(); iter != themap.end(); iter++) {
719 if (maxEnergy < (*iter).second) {
720 maxEnergy = (*iter).second;
721 unitWithMaxEnergy = (*iter).first;
725 return unitWithMaxEnergy;
729 int CentralEta,
int CentralPhi,
int CentralZ,
732 int goBackInEta = nCellInEta/2;
733 int goBackInPhi = nCellInPhi/2;
734 int matrixSize = nCellInEta*nCellInPhi;
735 crystalMatrix.clear();
736 crystalMatrix.resize(matrixSize);
738 int startEta = CentralZ*CentralEta - goBackInEta;
739 int startPhi = CentralPhi - goBackInPhi;
742 for (
int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
743 for(
int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
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(); }
749 crystalMatrix[i++] =
index;
756 int CentralX,
int CentralY,
int CentralZ,
758 int goBackInX = nCellInX/2;
759 int goBackInY = nCellInY/2;
760 crystalMatrix.clear();
762 int startX = CentralX - goBackInX;
763 int startY = CentralY - goBackInY;
765 for (
int ix = startX; ix < startX+nCellInX; ix ++ ) {
767 for(
int iy = startY; iy < startY + nCellInY; iy++ ) {
775 crystalMatrix.push_back(index);
T getParameter(std::string const &) const
EventNumber_t event() const
T getUntrackedParameter(std::string const &, T const &) const
~EcalRecHitsValidation()
Destructor.
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")
void analyze(const edm::Event &e, const edm::EventSetup &c)
Analyze.
std::vector< T >::const_iterator const_iterator
uint32_t getUnitWithMaxEnergy(MapType &themap)
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
int iphi() const
get the crystal iphi
uint32_t rawId() const
get the raw id
void findEndcapMatrix(int nCellInX, int nCellInY, int CentralX, int CentralY, int CentralZ, MapType &themap)
uint16_t getStatusCode() const
Tan< T >::type tan(const T &t)
EcalTrigTowerDetId tower() const
get the HCAL/trigger iphi of this crystal
int ieta() const
get the crystal ieta
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
const_iterator end() const
Log< T >::type log(const T &t)
std::vector< Item >::const_iterator const_iterator
T const * product() const
T const * product() const
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
int ietaAbs() const
get the absolute value of the crystal ieta
const_iterator begin() const
int zside() const
get the z-side of the crystal (1/-1)