38 EBuncalibrechitCollection_Token_ = consumes<EBUncalibratedRecHitCollection>(ps.
getParameter<
edm::InputTag>(
"EBuncalibrechitCollection"));
39 EEuncalibrechitCollection_Token_ = consumes<EEUncalibratedRecHitCollection>(ps.
getParameter<
edm::InputTag>(
"EEuncalibrechitCollection"));
48 if ( outputFile_.size() != 0 ) {
49 LogInfo(
"OutputInfo") <<
" Ecal RecHits Task histograms will be saved to '" << outputFile_.c_str() <<
"'";
51 LogInfo(
"OutputInfo") <<
" Ecal RecHits Task histograms will NOT be saved";
70 if ( verbose_ ) dbe_->showDirStructure();
78 meEBRecHitSimHitRatio_ = 0;
79 meEERecHitSimHitRatio_ = 0;
80 meESRecHitSimHitRatio_ = 0;
81 meEBRecHitSimHitRatio1011_ = 0;
82 meEERecHitSimHitRatio1011_ = 0;
83 meEBRecHitSimHitRatio12_ = 0;
84 meEERecHitSimHitRatio12_ = 0;
85 meEBRecHitSimHitRatio13_ = 0;
86 meEERecHitSimHitRatio13_ = 0;
87 meEBRecHitSimHitRatioGt35_ = 0;
88 meEERecHitSimHitRatioGt35_ = 0;
89 meEBUnRecHitSimHitRatio_ = 0;
90 meEEUnRecHitSimHitRatio_ = 0;
91 meEBUnRecHitSimHitRatioGt35_ = 0;
92 meEEUnRecHitSimHitRatioGt35_ = 0;
94 meEBe5x5OverSimHits_ = 0;
97 meEEe5x5OverSimHits_ = 0;
100 meEBRecHitLog10Energy_ = 0;
101 meEERecHitLog10Energy_ = 0;
102 meESRecHitLog10Energy_ = 0;
103 meEBRecHitLog10EnergyContr_ = 0;
104 meEERecHitLog10EnergyContr_ = 0;
105 meESRecHitLog10EnergyContr_ = 0;
106 meEBRecHitLog10Energy5x5Contr_ = 0;
107 meEERecHitLog10Energy5x5Contr_ = 0;
109 meEBRecHitsOccupancyFlag5_6_ = 0;
110 meEBRecHitsOccupancyFlag8_9_ = 0;
111 meEERecHitsOccupancyPlusFlag5_6_ = 0;
112 meEERecHitsOccupancyMinusFlag5_6_ = 0;
113 meEERecHitsOccupancyPlusFlag8_9_ = 0;
114 meEERecHitsOccupancyMinusFlag8_9_ = 0;
116 meEBRecHitFlags_ = 0;
117 meEBRecHitSimHitvsSimHitFlag5_6_ = 0;
118 meEBRecHitSimHitFlag6_ = 0;
119 meEBRecHitSimHitFlag7_ = 0;
120 meEB5x5RecHitSimHitvsSimHitFlag8_ = 0;
122 meEERecHitFlags_ = 0;
123 meEERecHitSimHitvsSimHitFlag5_6_ = 0;
124 meEERecHitSimHitFlag6_ = 0;
125 meEERecHitSimHitFlag7_ = 0;
132 dbe_->setCurrentFolder(
"EcalRecHitsV/EcalRecHitsTask");
134 histo =
"EcalRecHitsTask Gun Momentum";
135 meGunEnergy_ = dbe_->book1D(histo.c_str(), histo.c_str(), 100, 0., 1000.);
137 histo =
"EcalRecHitsTask Gun Eta";
138 meGunEta_ = dbe_->book1D(histo.c_str(), histo.c_str(), 700, -3.5, 3.5);
140 histo =
"EcalRecHitsTask Gun Phi";
141 meGunPhi_ = dbe_->book1D(histo.c_str(), histo.c_str(), 360, 0., 360.);
143 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio";
144 meEBRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
146 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio";
147 meEERecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
149 histo =
"EcalRecHitsTask Preshower RecSimHit Ratio";
150 meESRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
152 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio gt 3p5 GeV";
153 meEBRecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
155 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio gt 3p5 GeV";
156 meEERecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
158 histo =
"EcalRecHitsTask Barrel Unc RecSimHit Ratio";
159 meEBUnRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
161 histo =
"EcalRecHitsTask Endcap Unc RecSimHit Ratio";
162 meEEUnRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
164 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=10 11";
165 meEBRecHitSimHitRatio1011_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
167 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=10 11";
168 meEERecHitSimHitRatio1011_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
170 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=12";
171 meEBRecHitSimHitRatio12_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
173 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=12";
174 meEERecHitSimHitRatio12_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
176 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=13";
177 meEBRecHitSimHitRatio13_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
179 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=13";
180 meEERecHitSimHitRatio13_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
182 histo =
"EcalRecHitsTask Barrel Unc RecSimHit Ratio gt 3p5 GeV";
183 meEBUnRecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
185 histo =
"EcalRecHitsTask Endcap Unc RecSimHit Ratio gt 3p5 GeV";
186 meEEUnRecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
188 histo =
"EcalRecHitsTask Barrel Rec E5x5";
189 meEBe5x5_ = dbe_->book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
191 histo =
"EcalRecHitsTask Barrel Rec E5x5 over Sim E5x5";
192 meEBe5x5OverSimHits_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
194 histo =
"EcalRecHitsTask Barrel Rec E5x5 over gun energy";
195 meEBe5x5OverGun_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
197 histo =
"EcalRecHitsTask Endcap Rec E5x5";
198 meEEe5x5_ = dbe_->book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
200 histo =
"EcalRecHitsTask Endcap Rec E5x5 over Sim E5x5";
201 meEEe5x5OverSimHits_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
203 histo =
"EcalRecHitsTask Endcap Rec E5x5 over gun energy";
204 meEEe5x5OverGun_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
206 meEBRecHitLog10Energy_ = dbe_->book1D(
"EcalRecHitsTask Barrel Log10 Energy",
"EcalRecHitsTask Barrel Log10 Energy", 90, -5., 4. );
207 meEERecHitLog10Energy_ = dbe_->book1D(
"EcalRecHitsTask Endcap Log10 Energy",
"EcalRecHitsTask Endcap Log10 Energy", 90, -5., 4. );
208 meESRecHitLog10Energy_ = dbe_->book1D(
"EcalRecHitsTask Preshower Log10 Energy",
"EcalRecHitsTask Preshower Log10 Energy", 90, -5., 4. );
209 meEBRecHitLog10EnergyContr_ = dbe_->bookProfile(
"EcalRecHits Barrel Log10En vs Hit Contribution",
"EcalRecHits Barrel Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
210 meEERecHitLog10EnergyContr_ = dbe_->bookProfile(
"EcalRecHits Endcap Log10En vs Hit Contribution",
"EcalRecHits Endcap Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
211 meESRecHitLog10EnergyContr_ = dbe_->bookProfile(
"EcalRecHits Preshower Log10En vs Hit Contribution",
"EcalRecHits Preshower Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
212 meEBRecHitLog10Energy5x5Contr_ = dbe_->bookProfile(
"EcalRecHits Barrel Log10En5x5 vs Hit Contribution",
"EcalRecHits Barrel Log10En5x5 vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
213 meEERecHitLog10Energy5x5Contr_ = dbe_->bookProfile(
"EcalRecHits Endcap Log10En5x5 vs Hit Contribution",
"EcalRecHits Endcap Log10En5x5 vs Hit Contribution", 90, -5., 4., 100, 0., 1. );
216 histo =
"EB Occupancy Flag=5 6";
217 meEBRecHitsOccupancyFlag5_6_ = dbe_->book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
218 histo =
"EB Occupancy Flag=8 9";
219 meEBRecHitsOccupancyFlag8_9_ = dbe_->book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
221 histo =
"EE+ Occupancy Flag=5 6";
222 meEERecHitsOccupancyPlusFlag5_6_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
223 histo =
"EE- Occupancy Flag=5 6";
224 meEERecHitsOccupancyMinusFlag5_6_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
225 histo =
"EE+ Occupancy Flag=8 9";
226 meEERecHitsOccupancyPlusFlag8_9_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
227 histo =
"EE- Occupancy Flag=8 9";
228 meEERecHitsOccupancyMinusFlag8_9_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
231 histo =
"EcalRecHitsTask Barrel Reco Flags";
232 meEBRecHitFlags_ = dbe_->book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);
233 histo =
"EcalRecHitsTask Endcap Reco Flags";
234 meEERecHitFlags_ = dbe_->book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);
235 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio vs SimHit Flag=5 6";
236 meEBRecHitSimHitvsSimHitFlag5_6_ = dbe_->book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );
237 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio vs SimHit Flag=5 6";
238 meEERecHitSimHitvsSimHitFlag5_6_ = dbe_->book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );
239 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Flag=6";
240 meEBRecHitSimHitFlag6_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
241 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Flag=6";
242 meEERecHitSimHitFlag6_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
243 histo =
"EcalRecHitsTask Barrel RecSimHit Ratio Flag=7";
244 meEBRecHitSimHitFlag7_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
245 histo =
"EcalRecHitsTask Endcap RecSimHit Ratio Flag=7";
246 meEERecHitSimHitFlag7_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
247 histo =
"EcalRecHitsTask Barrel 5x5 RecSimHit Ratio vs SimHit Flag=8";
248 meEB5x5RecHitSimHitvsSimHitFlag8_ = dbe_->book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );
255 if ( outputFile_.size() != 0 &&
dbe_ )
dbe_->save(outputFile_);
272 LogInfo(
"EcalRecHitsTask, EventInfo: ") <<
" Run = " << e.
id().
run() <<
" Event = " << e.
id().
event();
278 const double barrelADCtoGeV_ = agc->
getEBValue();
279 const double endcapADCtoGeV_ = agc->
getEEValue();
284 if (!MCEvt.
isValid()) { skipMC =
true; }
288 bool skipBarrel =
false;
291 e.
getByToken( EBuncalibrechitCollection_Token_, EcalUncalibRecHitEB);
292 if (EcalUncalibRecHitEB.
isValid()) {
293 EBUncalibRecHit = EcalUncalibRecHitEB.
product() ;
298 bool skipEndcap =
false;
301 e.
getByToken( EEuncalibrechitCollection_Token_, EcalUncalibRecHitEE);
302 if (EcalUncalibRecHitEE.
isValid()){
303 EEUncalibRecHit = EcalUncalibRecHitEE.
product () ;
310 e.
getByToken( EBrechitCollection_Token_, EcalRecHitEB);
312 EBRecHit = EcalRecHitEB.
product();
319 e.
getByToken( EErechitCollection_Token_, EcalRecHitEE);
321 EERecHit = EcalRecHitEE.
product ();
326 bool skipPreshower =
false;
329 e.
getByToken( ESrechitCollection_Token_, EcalRecHitES);
331 ESRecHit = EcalRecHitES.
product ();
333 skipPreshower =
true;
341 for ( HepMC::GenEvent::particle_const_iterator
p = MCEvt->GetEvent()->particles_begin();
p != MCEvt->GetEvent()->particles_end(); ++
p ) {
342 double htheta = (*p)->momentum().theta();
343 double heta = -99999.;
344 if(
tan(htheta * 0.5) > 0 ) {
345 heta = -
log(
tan(htheta * 0.5));
347 double hphi = (*p)->momentum().phi();
348 hphi = (hphi>=0) ? hphi : hphi+2*
M_PI;
349 hphi = hphi /
M_PI * 180.;
351 LogDebug(
"EventInfo") <<
"EcalRecHitsTask: Particle gun type form MC = " <<
abs((*p)->pdg_id())
352 <<
"\n" <<
"Energy = "<< (*p)->momentum().e()
353 <<
"\n" <<
"Eta = " << heta
354 <<
"\n" <<
"Phi = " << hphi;
356 if ( (*p)->momentum().e() > eGun ) eGun = (*p)->momentum().e();
358 if (meGunEnergy_) meGunEnergy_->Fill((*p)->momentum().e());
359 if (meGunEta_) meGunEta_ ->Fill(heta);
360 if (meGunPhi_) meGunPhi_ ->Fill(hphi);
371 std::auto_ptr<MixCollection<PCaloHit> >
376 const int ebcSize = 90;
377 double ebcontr[ebcSize];
378 double ebcontr25[ebcSize];
379 for(
int i=0;
i<ebcSize;
i++ ) { ebcontr[
i] = 0.0; ebcontr25[
i] = 0.0; }
386 <<
"CaloHit " << hitItr->getName() <<
" DetID = " << hitItr->id() <<
"\n"
387 <<
"Energy = " << hitItr->energy() <<
" Time = " << hitItr->time() <<
"\n"
388 <<
"EBDetId = " << ebid.
ieta() <<
" " << ebid.
iphi();
390 uint32_t crystid = ebid.
rawId();
391 ebSimMap[crystid] += hitItr->energy();
402 if( myRecHit == EBRecHit->
end() )
continue;
403 ebRecMap[EBid.
rawId()] += myRecHit->energy();
406 ebtotal += myRecHit->energy();
407 if( myRecHit->energy() > 0 ) {
408 if( meEBRecHitLog10Energy_ ) meEBRecHitLog10Energy_->Fill( log10( myRecHit->energy() ) );
409 int log10i = int( ( log10( myRecHit->energy() ) + 5. ) * 10. );
410 if( log10i >=0 and log10i < ebcSize ) ebcontr[ log10i ] += myRecHit->energy();
414 if ( ebSimMap[EBid.
rawId()] != 0. ) {
415 double uncEnergy = uncalibRecHit->amplitude()*barrelADCtoGeV_;
416 if (meEBUnRecHitSimHitRatio_) {meEBUnRecHitSimHitRatio_ ->Fill(uncEnergy/ebSimMap[EBid.
rawId()]);}
417 if (meEBUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBUnRecHitSimHitRatioGt35_->Fill(uncEnergy/ebSimMap[EBid.
rawId()]);}
420 if (myRecHit != EBRecHit->
end()) {
421 if ( ebSimMap[EBid.
rawId()] != 0. ) {
422 if (meEBRecHitSimHitRatio_) {meEBRecHitSimHitRatio_ ->Fill(myRecHit->energy()/ebSimMap[EBid.
rawId()]);}
423 if (meEBRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBRecHitSimHitRatioGt35_->Fill(myRecHit->energy()/ebSimMap[EBid.
rawId()]);}
432 if( csmi != ecs->
end() ) csc = *csmi;
436 if( meEBRecHitSimHitRatio1011_ != 0 &&
437 ( sc == 10 || sc == 11 ) ) { meEBRecHitSimHitRatio1011_->Fill(myRecHit->energy()/ebSimMap[EBid.
rawId()]); }
438 if( meEBRecHitSimHitRatio12_ != 0 && sc == 12 ) { meEBRecHitSimHitRatio12_->Fill(myRecHit->energy()/ebSimMap[EBid.
rawId()]); }
444 double ttSimEnergy = 0;
448 for( std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); dit++ ) {
450 ttSimEnergy += ebSimMap[ttEBid.
rawId()];
452 if( vid.size() != 0 ) ttSimEnergy = ttSimEnergy / vid.size();
454 if( meEBRecHitSimHitRatio13_ != 0 && sc == 13 && ttSimEnergy != 0 )
455 meEBRecHitSimHitRatio13_->Fill(myRecHit->energy()/ttSimEnergy);
457 int flag = myRecHit->recoFlag();
458 if( meEBRecHitFlags_ != 0 ) meEBRecHitFlags_->Fill( flag );
460 meEBRecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/ebSimMap[EBid.
rawId()], ebSimMap[EBid.
rawId()] );
462 meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.
rawId()] );
464 meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.
rawId()] );
466 meEB5x5RecHitSimHitvsSimHitFlag8_->Fill( myRecHit->energy()/ttSimEnergy, ttSimEnergy );
469 meEBRecHitsOccupancyFlag5_6_ ->
Fill(EBid.
ieta(), EBid.
iphi());
471 meEBRecHitsOccupancyFlag8_9_ ->
Fill(EBid.
ieta(), EBid.
iphi());
480 uint32_t ebcenterid = getUnitWithMaxEnergy(ebRecMap);
483 int by = myEBid.
iphi();
484 int bz = myEBid.
zside();
485 findBarrelMatrix(5,5,bx,by,bz,ebRecMap);
488 for (
unsigned int i = 0;
i < crystalMatrix.size();
i++ ) {
489 e5x5rec += ebRecMap[crystalMatrix[
i]];
490 e5x5sim += ebSimMap[crystalMatrix[
i]];
491 if( ebRecMap[crystalMatrix[
i]] > 0 ) {
492 int log10i25 = int( ( log10( ebRecMap[crystalMatrix[
i]] ) + 5. ) * 10. );
493 if( log10i25 >=0 && log10i25 < ebcSize ) ebcontr25[ log10i25 ] += ebRecMap[crystalMatrix[
i]];
497 if( meEBe5x5_ ) meEBe5x5_->Fill(e5x5rec);
498 if ( e5x5sim > 0. && meEBe5x5OverSimHits_ ) meEBe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
499 if ( eGun > 0. && meEBe5x5OverGun_ ) meEBe5x5OverGun_->Fill(e5x5rec/eGun);
501 if( meEBRecHitLog10EnergyContr_ && ebtotal != 0 ) {
502 for(
int i=0;
i<ebcSize;
i++ ) {
503 meEBRecHitLog10EnergyContr_->Fill( -5.+(
float(
i)+0.5)/10., ebcontr[
i]/ebtotal );
507 if( meEBRecHitLog10Energy5x5Contr_ && e5x5rec != 0 ) {
508 for(
int i=0;
i<ebcSize;
i++ ) {
509 meEBRecHitLog10Energy5x5Contr_->Fill( -5.+(
float(
i)+0.5)/10., ebcontr25[
i]/e5x5rec );
517 if ( ! skipEndcap ) {
521 std::auto_ptr<MixCollection<PCaloHit> >
526 const int eecSize = 90;
527 double eecontr[eecSize];
528 double eecontr25[eecSize];
529 for(
int i=0;
i<eecSize;
i++ ) { eecontr[
i] = 0.0; eecontr25[
i] = 0.0; }
536 <<
" CaloHit " << hitItr->getName() <<
" DetID = " << hitItr->id() <<
"\n"
537 <<
"Energy = " << hitItr->energy() <<
" Time = " << hitItr->time() <<
"\n"
538 <<
"EEDetId side " << eeid.
zside() <<
" = " << eeid.
ix() <<
" " << eeid.
iy();
540 uint32_t crystid = eeid.
rawId();
541 eeSimMap[crystid] += hitItr->energy();
552 if( myRecHit == EERecHit->
end() )
continue;
553 eeRecMap[EEid.
rawId()] += myRecHit->energy();
556 eetotal += myRecHit->energy();
557 if( myRecHit->energy() > 0 ) {
558 if( meEERecHitLog10Energy_ ) meEERecHitLog10Energy_->Fill( log10( myRecHit->energy() ) );
559 int log10i = int( ( log10( myRecHit->energy() ) + 5. ) * 10. );
560 if( log10i >=0 and log10i < eecSize ) eecontr[ log10i ] += myRecHit->energy();
564 if ( eeSimMap[EEid.
rawId()] != 0. ) {
565 double uncEnergy = uncalibRecHit->amplitude()*endcapADCtoGeV_;
566 if (meEEUnRecHitSimHitRatio_) {meEEUnRecHitSimHitRatio_ ->Fill(uncEnergy/eeSimMap[EEid.
rawId()]);}
567 if (meEEUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEEUnRecHitSimHitRatioGt35_->Fill(uncEnergy/eeSimMap[EEid.
rawId()]);}
570 if (myRecHit != EERecHit->
end()) {
571 if ( eeSimMap[EEid.
rawId()] != 0. ) {
572 if (meEERecHitSimHitRatio_) {meEERecHitSimHitRatio_ ->Fill(myRecHit->energy()/eeSimMap[EEid.
rawId()]); }
573 if (meEERecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEERecHitSimHitRatioGt35_->Fill(myRecHit->energy()/eeSimMap[EEid.
rawId()]); }
582 if( csmi != ecs->
end() ) csc = *csmi;
584 if( meEERecHitSimHitRatio1011_ != 0 &&
585 ( sc == 10 || sc == 11 ) ) { meEERecHitSimHitRatio1011_->Fill(myRecHit->energy()/eeSimMap[EEid.
rawId()]); }
586 if( meEERecHitSimHitRatio12_ != 0 && sc == 12 ) { meEERecHitSimHitRatio12_->Fill(myRecHit->energy()/eeSimMap[EEid.
rawId()]); }
587 if( meEERecHitSimHitRatio13_ != 0 && sc == 13 ) { meEERecHitSimHitRatio13_->Fill(myRecHit->energy()/eeSimMap[EEid.
rawId()]); }
590 int flag = myRecHit->recoFlag();
591 if( meEERecHitFlags_ != 0 ) meEERecHitFlags_->Fill( flag );
593 meEERecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/eeSimMap[EEid.
rawId()], eeSimMap[EEid.
rawId()] );
595 meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.
rawId()] );
597 meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.
rawId()] );
599 if (EEid.
zside() > 0) {
601 meEERecHitsOccupancyPlusFlag5_6_ ->
Fill(EEid.
ix(), EEid.
iy());
603 meEERecHitsOccupancyPlusFlag8_9_ ->
Fill(EEid.
ix(), EEid.
iy());
605 if (EEid.
zside() < 0) {
607 meEERecHitsOccupancyMinusFlag5_6_ ->
Fill(EEid.
ix(), EEid.
iy());
609 meEERecHitsOccupancyMinusFlag8_9_ ->
Fill(EEid.
ix(), EEid.
iy());
618 uint32_t eecenterid = getUnitWithMaxEnergy(eeRecMap);
620 int bx = myEEid.
ix();
621 int by = myEEid.
iy();
622 int bz = myEEid.
zside();
623 findEndcapMatrix(5,5,bx,by,bz,eeRecMap);
626 for (
unsigned int i = 0;
i < crystalMatrix.size();
i++ ) {
627 e5x5rec += eeRecMap[crystalMatrix[
i]];
628 e5x5sim += eeSimMap[crystalMatrix[
i]];
629 if( eeRecMap[crystalMatrix[
i]] > 0 ) {
630 int log10i25 = int( ( log10( eeRecMap[crystalMatrix[
i]] ) + 5. ) * 10. );
631 if( log10i25 >=0 && log10i25 < eecSize ) eecontr25[ log10i25 ] += eeRecMap[crystalMatrix[
i]];
635 if( meEEe5x5_ ) meEEe5x5_->Fill(e5x5rec);
636 if ( e5x5sim > 0. && meEEe5x5OverSimHits_ ) meEEe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
637 if ( eGun > 0. && meEEe5x5OverGun_ ) meEEe5x5OverGun_->Fill(e5x5rec/eGun);
640 if( meEERecHitLog10EnergyContr_ && eetotal != 0 ) {
641 for(
int i=0;
i<eecSize;
i++ ) {
642 meEERecHitLog10EnergyContr_->Fill( -5.+(
float(
i)+0.5)/10., eecontr[
i]/eetotal );
646 if( meEERecHitLog10Energy5x5Contr_ && e5x5rec != 0 ) {
647 for(
int i=0;
i<eecSize;
i++ ) {
648 meEERecHitLog10Energy5x5Contr_->Fill( -5.+(
float(
i)+0.5)/10., eecontr25[
i]/e5x5rec );
656 if ( ! skipPreshower ) {
660 std::auto_ptr<MixCollection<PCaloHit> >
664 const int escSize = 90;
665 double escontr[escSize];
666 for(
int i=0;
i<escSize;
i++ ) { escontr[
i] = 0.0; }
674 <<
" CaloHit " << hitItr->getName() <<
" DetID = " << hitItr->id() <<
"\n"
675 <<
"Energy = " << hitItr->energy() <<
" Time = " << hitItr->time() <<
"\n"
676 <<
"ESDetId strip " << esid.
strip() <<
" = " << esid.
six() <<
" " << esid.
siy();
678 uint32_t crystid = esid.
rawId();
679 esSimMap[crystid] += hitItr->energy();
686 if ( esSimMap[ESid.
rawId()] != 0. ) {
689 estotal += recHit->energy();
690 if( recHit->energy() > 0 ) {
691 if( meESRecHitLog10Energy_ ) meESRecHitLog10Energy_->Fill( log10( recHit->energy() ) );
692 int log10i = int( ( log10( recHit->energy() ) + 5. ) * 10. );
693 if( log10i >=0 and log10i < escSize ) escontr[ log10i ] += recHit->energy();
696 if (meESRecHitSimHitRatio_) {
697 meESRecHitSimHitRatio_ ->Fill(recHit->energy()/esSimMap[ESid.
rawId()]);
704 if( meESRecHitLog10EnergyContr_ && estotal != 0 ) {
705 for(
int i=0;
i<escSize;
i++ ) {
706 meESRecHitLog10EnergyContr_->Fill( -5.+(
float(
i)+0.5)/10., escontr[
i]/estotal );
719 uint32_t unitWithMaxEnergy = 0;
720 float maxEnergy = 0.;
722 MapType::iterator
iter;
723 for (iter = themap.begin(); iter != themap.end(); iter++) {
725 if (maxEnergy < (*iter).second) {
726 maxEnergy = (*iter).second;
727 unitWithMaxEnergy = (*iter).first;
731 return unitWithMaxEnergy;
735 int CentralEta,
int CentralPhi,
int CentralZ,
738 int goBackInEta = nCellInEta/2;
739 int goBackInPhi = nCellInPhi/2;
740 int matrixSize = nCellInEta*nCellInPhi;
741 crystalMatrix.clear();
742 crystalMatrix.resize(matrixSize);
744 int startEta = CentralZ*CentralEta - goBackInEta;
745 int startPhi = CentralPhi - goBackInPhi;
748 for (
int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
749 for(
int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
751 if (
abs(ieta) > 85 ||
abs(ieta)<1 ) {
continue; }
752 if (iphi< 1) { index =
EBDetId(ieta,iphi+360).
rawId(); }
753 else if(iphi>360) { index =
EBDetId(ieta,iphi-360).
rawId(); }
755 crystalMatrix[i++] =
index;
762 int CentralX,
int CentralY,
int CentralZ,
764 int goBackInX = nCellInX/2;
765 int goBackInY = nCellInY/2;
766 crystalMatrix.clear();
768 int startX = CentralX - goBackInX;
769 int startY = CentralY - goBackInY;
771 for (
int ix = startX; ix < startX+nCellInX; ix ++ ) {
773 for(
int iy = startY; iy < startY + nCellInY; iy++ ) {
781 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.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
void analyze(const edm::Event &e, const edm::EventSetup &c)
Analyze.
Code getStatusCode() const
return decoded status
std::vector< EcalUncalibratedRecHit >::const_iterator const_iterator
uint32_t getUnitWithMaxEnergy(MapType &themap)
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)
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Tan< T >::type tan(const T &t)
Abs< T >::type abs(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.
const_iterator end() const
T const * product() const
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
std::vector< Item >::const_iterator const_iterator
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)