18 g4InfoLabel(ps.getParameter<std::
string>(
"moduleLabelG4")),
19 EEHitsCollection(ps.getParameter<std::
string>(
"EEHitsCollection")),
20 ValidationCollection(ps.getParameter<std::
string>(
"ValidationCollection")){
29 if (
verbose_ ) { dbe_->setVerbose(1); }
30 else { dbe_->setVerbose(0); }
34 if (
verbose_ ) dbe_->showDirStructure();
68 for (
int myStep = 0; myStep<26; myStep++) {
eRLength[myStep] = 0.0; }
73 dbe_->setCurrentFolder(
"EcalHitsV/EcalSimHitsValidation");
75 sprintf (histo,
"EE+ hits multiplicity" ) ;
76 meEEzpHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
78 sprintf (histo,
"EE- hits multiplicity" ) ;
79 meEEzmHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
81 sprintf (histo,
"EE+ crystals multiplicity" ) ;
84 sprintf (histo,
"EE- crystals multiplicity" ) ;
87 sprintf (histo,
"EE+ occupancy" ) ;
90 sprintf (histo,
"EE- occupancy" ) ;
93 sprintf (histo,
"EE longitudinal shower profile" ) ;
96 sprintf (histo,
"EE hits energy spectrum" );
99 sprintf (histo,
"EE hits log10energy spectrum" );
102 sprintf (histo,
"EE hits log10energy spectrum vs normalized energy" );
105 sprintf (histo,
"EE hits log10energy spectrum vs normalized energy25" );
108 sprintf (histo,
"EE hits energy spectrum 2" );
111 sprintf (histo,
"EE crystal energy spectrum" );
114 sprintf (histo,
"EE crystal energy spectrum 2" );
117 sprintf (histo,
"EE E1" ) ;
118 meEEe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
120 sprintf (histo,
"EE E4" ) ;
121 meEEe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
123 sprintf (histo,
"EE E9" ) ;
124 meEEe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
126 sprintf (histo,
"EE E16" ) ;
127 meEEe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
129 sprintf (histo,
"EE E25" ) ;
130 meEEe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
132 sprintf (histo,
"EE E1oE4" ) ;
133 meEEe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
135 sprintf (histo,
"EE E1oE9" ) ;
136 meEEe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
138 sprintf (histo,
"EE E4oE9" ) ;
139 meEEe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
141 sprintf (histo,
"EE E9oE16" ) ;
142 meEEe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
144 sprintf (histo,
"EE E1oE25" ) ;
145 meEEe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
147 sprintf (histo,
"EE E9oE25" ) ;
148 meEEe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
150 sprintf (histo,
"EE E16oE25" ) ;
151 meEEe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
179 if( ! EcalHitsEE.
isValid() )
return;
184 std::vector<PCaloHit> theEECaloHits;
185 theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
189 std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
199 std::vector<double> econtr(140, 0. );
200 std::vector<double> econtr25(140, 0. );
205 uint32_t nEEzpHits = 0;
206 uint32_t nEEzmHits = 0;
208 for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin();
209 isim != theEECaloHits.end(); ++isim){
211 if ( isim->time() > 500. ) {
continue; }
213 CaloHitMap[ isim->id()].push_back(&(*isim));
218 <<
" CaloHit " << isim->getName() <<
"\n"
219 <<
" DetID = "<<isim->id()<<
" EEDetId = " << eeid.ix() <<
" " << eeid.iy() <<
"\n"
220 <<
" Time = " << isim->time() <<
"\n"
221 <<
" Track Id = " << isim->geantTrackId() <<
"\n"
222 <<
" Energy = " << isim->energy();
224 uint32_t crystid = eeid.rawId();
226 if (eeid.zside() > 0 ) {
228 EEetzp_ += isim->energy();
229 eemapzp[crystid] += isim->energy();
232 else if (eeid.zside() < 0 ) {
234 EEetzm_ += isim->energy();
235 eemapzm[crystid] += isim->energy();
240 if( isim->energy() > 0 ) {
242 int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
243 if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
246 eemap[crystid] += isim->energy();
255 for (
std::map<uint32_t,
float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it )
meEEcrystalEnergy_->
Fill((*it).second);
265 int nEEHits = nEEzmHits + nEEzpHits;
270 int bx = myEEid.
ix();
271 int by = myEEid.
iy();
272 int bz = myEEid.
zside();
280 std::vector<uint32_t> ids25; ids25 =
getIdsAroundMax(5,5,bx,by,bz,eemap);
282 for(
unsigned i=0;
i<25;
i++ ) {
283 for(
unsigned int j=0;
j<CaloHitMap[ids25[
i]].size();
j++ ) {
284 if( CaloHitMap[ids25[
i]][
j]->
energy() > 0 ) {
285 int log10i = int( ( log10( CaloHitMap[ids25[
i]][
j]->
energy()) + 10. ) * 10. );
286 if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[
i]][
j]->energy();
310 for(
int i=0;
i<140;
i++ ) {
316 for(
int i=0;
i<140;
i++ ) {
323 if( MyPEcalValidInfo.
isValid() ) {
324 if ( MyPEcalValidInfo->ee1x1() > 0. ) {
325 std::vector<float> EX0 = MyPEcalValidInfo->eX0();
327 for (
int myStep=0; myStep< 26; myStep++ ) {
336 int centralX,
int centralY,
337 int centralZ,
MapType& themap){
340 float totalEnergy = 0.;
342 int goBackInX = nCellInX/2;
343 int goBackInY = nCellInY/2;
344 int startX = centralX-goBackInX;
345 int startY = centralY-goBackInY;
347 for (
int ix=startX; ix<startX+nCellInX; ix++) {
348 for (
int iy=startY; iy<startY+nCellInY; iy++) {
356 totalEnergy += themap[
index];
362 << nCellInX <<
" x " << nCellInY
363 <<
" EE matrix energy = " << totalEnergy
364 <<
" for " << ncristals <<
" crystals";
372 std::vector<uint32_t> ids( nCellInX*nCellInY );
374 int goBackInX = nCellInX/2;
375 int goBackInY = nCellInY/2;
376 int startX = centralX-goBackInX;
377 int startY = centralY-goBackInY;
379 for (
int ix=startX; ix<startX+nCellInX; ix++) {
380 for (
int iy=startY; iy<startY+nCellInY; iy++) {
388 ids[ncristals] =
index;
397 int CentralX,
int CentralY,
int CentralZ,
400 int goBackInX = nCellInX/2;
401 int goBackInY = nCellInY/2;
403 int startX = CentralX - goBackInX;
404 int startY = CentralY - goBackInY;
407 for (
int ix = startX; ix < startX+nCellInX; ix ++ ) {
409 for(
int iy = startY; iy < startY + nCellInY; iy++ ) {
417 fillmap[i++] = themap[
index];
422 if ( fillmap[i/2] == themap[centerid] )
431 float e012 = themap[0]+themap[1]+themap[2];
432 float e036 = themap[0]+themap[3]+themap[6];
433 float e678 = themap[6]+themap[7]+themap[8];
434 float e258 = themap[2]+themap[5]+themap[8];
436 if ( (e012>e678 || e012==e678) && (e036>e258 || e036==e258))
437 return E22=themap[0]+themap[1]+themap[3]+themap[4];
438 else if ( (e012>e678 || e012==e678) && (e036<e258 || e036==e258) )
439 return E22=themap[1]+themap[2]+themap[4]+themap[5];
440 else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
441 return E22=themap[3]+themap[4]+themap[6]+themap[7];
442 else if ( (e012<e678|| e012==e678) && (e036<e258|| e036==e258) )
443 return E22=themap[4]+themap[5]+themap[7]+themap[8];
451 float e0_4 = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
452 float e0_20 = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
453 float e4_24 = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
454 float e0_24 = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
456 if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
457 return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
458 else if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
459 return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
460 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
461 return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
462 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
463 return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
472 uint32_t unitWithMaxEnergy = 0;
473 float maxEnergy = 0.;
475 MapType::iterator iter;
476 for (iter = themap.begin(); iter != themap.end(); iter++) {
478 if (maxEnergy < (*iter).second) {
479 maxEnergy = (*iter).second;
480 unitWithMaxEnergy = (*iter).first;
485 <<
" max energy of " << maxEnergy
486 <<
" GeV in Unit id " << unitWithMaxEnergy;
487 return unitWithMaxEnergy;
EventNumber_t event() const
MonitorElement * meEEe9oe16_
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * meEEzpCrystals_
virtual float energyInMatrixEE(int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap)
MonitorElement * meEEe16_
MonitorElement * meEEzmHits_
MonitorElement * meEEe1oe4_
MonitorElement * meEEzmOccupancy_
MonitorElement * meEEHitEnergy2_
~EcalEndcapSimHitsValidation()
Destructor.
uint32_t rawId() const
get the raw id
MonitorElement * meEEcrystalEnergy_
MonitorElement * meEEe4oe9_
MonitorElement * meEEzpOccupancy_
MonitorElement * meEEe16oe25_
float eCluster4x4(float e33, MapType &themap)
std::string EEHitsCollection
MonitorElement * meEELongitudinalShower_
std::string ValidationCollection
MonitorElement * meEEe9oe25_
MonitorElement * meEEzmCrystals_
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
MonitorElement * meEEhitLog10Energy25Norm_
bool fillEEMatrix(int nCellInX, int nCellInY, int CentralX, int CentralY, int CentralZ, MapType &fillmap, MapType &themap)
MonitorElement * meEEcrystalEnergy2_
MonitorElement * meEEHitEnergy_
MonitorElement * meEEzpHits_
MonitorElement * meEEe1oe9_
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
MonitorElement * meEEe1oe25_
void analyze(const edm::Event &e, const edm::EventSetup &c)
Analyze.
uint32_t getUnitWithMaxEnergy(MapType &themap)
std::map< uint32_t, float, std::less< uint32_t > > MapType
MonitorElement * meEEe25_
MonitorElement * meEEhitLog10Energy_
MonitorElement * meEEhitLog10EnergyNorm_
void Reset(void)
reset ME (ie. contents, errors, etc)
std::vector< uint32_t > getIdsAroundMax(int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap)
float eCluster2x2(MapType &themap)
EcalEndcapSimHitsValidation(const edm::ParameterSet &ps)
Constructor.