18 g4InfoLabel(ps.getParameter<std::
string>(
"moduleLabelG4")),
19 EBHitsCollection(ps.getParameter<std::
string>(
"EBHitsCollection")),
20 ValidationCollection(ps.getParameter<std::
string>(
"ValidationCollection")){
29 if (
verbose_ ) { dbe_->setVerbose(1); }
30 else { dbe_->setVerbose(0); }
34 if (
verbose_ ) dbe_->showDirStructure();
64 for (
int myStep = 0; myStep<26; myStep++) {
eRLength[myStep] = 0.0; }
69 dbe_->setCurrentFolder(
"EcalHitsV/EcalSimHitsValidation");
71 sprintf (histo,
"EB hits multiplicity" ) ;
72 menEBHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
74 sprintf (histo,
"EB crystals multiplicity" ) ;
77 sprintf (histo,
"EB occupancy" ) ;
78 meEBOccupancy_ = dbe_->book2D(histo, histo, 360, 0., 360., 170, -85., 85.);
80 sprintf (histo,
"EB longitudinal shower profile" ) ;
83 sprintf (histo,
"EB hits energy spectrum" );
86 sprintf (histo,
"EB hits log10energy spectrum" );
89 sprintf (histo,
"EB hits log10energy spectrum vs normalized energy" );
92 sprintf (histo,
"EB hits log10energy spectrum vs normalized energy25" );
95 sprintf (histo,
"EB hits energy spectrum 2" );
98 sprintf (histo,
"EB crystal energy spectrum" );
101 sprintf (histo,
"EB crystal energy spectrum 2" );
104 sprintf (histo,
"EB E1" ) ;
105 meEBe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
107 sprintf (histo,
"EB E4" ) ;
108 meEBe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
110 sprintf (histo,
"EB E9" ) ;
111 meEBe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
113 sprintf (histo,
"EB E16" ) ;
114 meEBe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
116 sprintf (histo,
"EB E25" ) ;
117 meEBe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
119 sprintf (histo,
"EB E1oE4" ) ;
120 meEBe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
122 sprintf (histo,
"EB E1oE9" ) ;
123 meEBe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
125 sprintf (histo,
"EB E4oE9" ) ;
126 meEBe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
128 sprintf (histo,
"EB E9oE16" ) ;
129 meEBe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
131 sprintf (histo,
"EB E1oE25" ) ;
132 meEBe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
134 sprintf (histo,
"EB E9oE25" ) ;
135 meEBe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
137 sprintf (histo,
"EB E16oE25" ) ;
138 meEBe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
165 if( ! EcalHitsEB.
isValid() )
return;
170 std::vector<PCaloHit> theEBCaloHits;
171 theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
175 double EBEnergy_ = 0.;
176 std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
183 std::vector<double> econtr(140, 0. );
184 std::vector<double> econtr25(140, 0. );
187 uint32_t nEBHits = 0;
189 for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin();
190 isim != theEBCaloHits.end(); ++isim){
192 if ( isim->time() > 500. ) {
continue; }
194 CaloHitMap[ isim->id()].push_back( &(*isim) );
199 <<
" CaloHit " << isim->getName() <<
"\n"
200 <<
" DetID = " << isim->id() <<
" EBDetId = " << ebid.ieta() <<
" " << ebid.iphi() <<
"\n"
201 <<
" Time = " << isim->time() <<
"\n"
202 <<
" Track Id = " << isim->geantTrackId() <<
"\n"
203 <<
" Energy = " << isim->energy();
207 uint32_t crystid = ebid.rawId();
208 ebmap[crystid] += isim->energy();
210 EBEnergy_ += isim->energy();
213 if( isim->energy() > 0 ) {
215 int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
216 if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
224 for (
std::map<uint32_t,
float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it )
meEBcrystalEnergy_->
Fill((*it).second);
237 int by = myEBid.
iphi();
238 int bz = myEBid.
zside();
246 std::vector<uint32_t> ids25; ids25 =
getIdsAroundMax(5,5,bx,by,bz,ebmap);
248 for(
unsigned i=0;
i<25;
i++ ) {
249 for(
unsigned int j=0;
j<CaloHitMap[ids25[
i]].size();
j++ ) {
250 if( CaloHitMap[ids25[
i]][
j]->
energy() > 0 ) {
251 int log10i = int( ( log10( CaloHitMap[ids25[
i]][
j]->
energy()) + 10. ) * 10. );
252 if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[
i]][
j]->energy();
276 for(
int i=0;
i<140;
i++ ) {
282 for(
int i=0;
i<140;
i++ ) {
291 if( MyPEcalValidInfo.
isValid() ) {
292 if ( MyPEcalValidInfo->eb1x1() > 0. ) {
293 std::vector<float> BX0 = MyPEcalValidInfo->bX0();
295 for (
int myStep=0; myStep< 26; myStep++ ) {
304 int centralEta,
int centralPhi,
305 int centralZ,
MapType& themap){
308 float totalEnergy = 0.;
310 int goBackInEta = nCellInEta/2;
311 int goBackInPhi = nCellInPhi/2;
312 int startEta = centralZ*centralEta-goBackInEta;
313 int startPhi = centralPhi-goBackInPhi;
315 for (
int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
316 for (
int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
319 if (
abs(ieta) > 85 ||
abs(ieta)<1 ) {
continue; }
320 if (iphi< 1) { index =
EBDetId(ieta,iphi+360).
rawId(); }
321 else if(iphi>360) { index =
EBDetId(ieta,iphi-360).
rawId(); }
324 totalEnergy += themap[
index];
330 << nCellInEta <<
" x " << nCellInPhi
331 <<
" EB matrix energy = " << totalEnergy
332 <<
" for " << ncristals <<
" crystals" ;
340 std::vector<uint32_t> ids( nCellInEta*nCellInPhi );
342 int goBackInEta = nCellInEta/2;
343 int goBackInPhi = nCellInPhi/2;
344 int startEta = centralZ*centralEta-goBackInEta;
345 int startPhi = centralPhi-goBackInPhi;
347 for (
int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
348 for (
int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
351 if (
abs(ieta) > 85 ||
abs(ieta)<1 ) {
continue; }
352 if (iphi< 1) { index =
EBDetId(ieta,iphi+360).
rawId(); }
353 else if(iphi>360) { index =
EBDetId(ieta,iphi-360).
rawId(); }
355 ids[ncristals] =
index;
365 int CentralEta,
int CentralPhi,
int CentralZ,
368 int goBackInEta = nCellInEta/2;
369 int goBackInPhi = nCellInPhi/2;
371 int startEta = CentralZ*CentralEta - goBackInEta;
372 int startPhi = CentralPhi - goBackInPhi;
375 for (
int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
376 for(
int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
378 if (
abs(ieta) > 85 ||
abs(ieta)<1 ) {
continue; }
379 if (iphi< 1) { index =
EBDetId(ieta,iphi+360).
rawId(); }
380 else if(iphi>360) { index =
EBDetId(ieta,iphi-360).
rawId(); }
382 fillmap[i++] = themap[
index];
388 if ( fillmap[i/2] == themap[ebcenterid] )
396 float e012 = themap[0]+themap[1]+themap[2];
397 float e036 = themap[0]+themap[3]+themap[6];
398 float e678 = themap[6]+themap[7]+themap[8];
399 float e258 = themap[2]+themap[5]+themap[8];
401 if ( (e012>e678 || e012==e678) && (e036>e258 || e036==e258))
402 return E22=themap[0]+themap[1]+themap[3]+themap[4];
403 else if ( (e012>e678 || e012==e678) && (e036<e258 || e036==e258) )
404 return E22=themap[1]+themap[2]+themap[4]+themap[5];
405 else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
406 return E22=themap[3]+themap[4]+themap[6]+themap[7];
407 else if ( (e012<e678|| e012==e678) && (e036<e258|| e036==e258) )
408 return E22=themap[4]+themap[5]+themap[7]+themap[8];
416 float e0_4 = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
417 float e0_20 = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
418 float e4_24 = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
419 float e0_24 = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
421 if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
422 return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
423 else if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
424 return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
425 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
426 return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
427 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
428 return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
437 uint32_t unitWithMaxEnergy = 0;
438 float maxEnergy = 0.;
440 MapType::iterator iter;
441 for (iter = themap.begin(); iter != themap.end(); iter++) {
443 if (maxEnergy < (*iter).second) {
444 maxEnergy = (*iter).second;
445 unitWithMaxEnergy = (*iter).first;
450 <<
" max energy of " << maxEnergy
451 <<
" GeV in Unit id " << unitWithMaxEnergy;
452 return unitWithMaxEnergy;
EventNumber_t event() const
T getUntrackedParameter(std::string const &, T const &) const
std::map< uint32_t, float, std::less< uint32_t > > MapType
std::vector< uint32_t > getIdsAroundMax(int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType &themap)
MonitorElement * meEBhitLog10Energy_
MonitorElement * meEBe16oe25_
MonitorElement * meEBe16_
void analyze(const edm::Event &e, const edm::EventSetup &c)
Analyze.
MonitorElement * meEBe9oe16_
MonitorElement * meEBe1oe9_
~EcalBarrelSimHitsValidation()
Destructor.
MonitorElement * meEBhitLog10EnergyNorm_
MonitorElement * meEBe9oe25_
MonitorElement * menEBCrystals_
MonitorElement * meEBe25_
MonitorElement * meEBhitLog10Energy25Norm_
MonitorElement * menEBHits_
int iphi() const
get the crystal iphi
uint32_t rawId() const
get the raw id
float eCluster4x4(float e33, MapType &themap)
MonitorElement * meEBe1oe4_
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
std::string ValidationCollection
MonitorElement * meEBLongitudinalShower_
Abs< T >::type abs(const T &t)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
MonitorElement * meEBe4oe9_
uint32_t getUnitWithMaxEnergy(MapType &themap)
virtual float energyInMatrixEB(int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType &themap)
MonitorElement * meEBcrystalEnergy2_
std::string EBHitsCollection
MonitorElement * meEBhitEnergy2_
float eCluster2x2(MapType &themap)
bool fillEBMatrix(int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &fillmap, MapType &themap)
MonitorElement * meEBOccupancy_
void Reset(void)
reset ME (ie. contents, errors, etc)
int ietaAbs() const
get the absolute value of the crystal ieta
MonitorElement * meEBcrystalEnergy_
MonitorElement * meEBhitEnergy_
int zside() const
get the z-side of the crystal (1/-1)
EcalBarrelSimHitsValidation(const edm::ParameterSet &ps)
Constructor.
MonitorElement * meEBe1oe25_