18 g4InfoLabel(ps.getParameter<std::
string>(
"moduleLabelG4")),
19 EEHitsCollection(ps.getParameter<std::
string>(
"EEHitsCollection")),
20 ValidationCollection(ps.getParameter<std::
string>(
"ValidationCollection")){
31 if ( verbose_ ) { dbe_->setVerbose(1); }
32 else { dbe_->setVerbose(0); }
36 if ( verbose_ ) dbe_->showDirStructure();
70 for (
int myStep = 0; myStep<26; myStep++) {
eRLength[myStep] = 0.0; }
75 dbe_->setCurrentFolder(
"EcalHitsV/EcalSimHitsValidation");
77 sprintf (histo,
"EE+ hits multiplicity" ) ;
78 meEEzpHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
80 sprintf (histo,
"EE- hits multiplicity" ) ;
81 meEEzmHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
83 sprintf (histo,
"EE+ crystals multiplicity" ) ;
86 sprintf (histo,
"EE- crystals multiplicity" ) ;
89 sprintf (histo,
"EE+ occupancy" ) ;
92 sprintf (histo,
"EE- occupancy" ) ;
95 sprintf (histo,
"EE longitudinal shower profile" ) ;
98 sprintf (histo,
"EE hits energy spectrum" );
101 sprintf (histo,
"EE hits log10energy spectrum" );
104 sprintf (histo,
"EE hits log10energy spectrum vs normalized energy" );
107 sprintf (histo,
"EE hits log10energy spectrum vs normalized energy25" );
110 sprintf (histo,
"EE hits energy spectrum 2" );
113 sprintf (histo,
"EE crystal energy spectrum" );
116 sprintf (histo,
"EE crystal energy spectrum 2" );
119 sprintf (histo,
"EE E1" ) ;
120 meEEe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
122 sprintf (histo,
"EE E4" ) ;
123 meEEe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
125 sprintf (histo,
"EE E9" ) ;
126 meEEe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
128 sprintf (histo,
"EE E16" ) ;
129 meEEe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
131 sprintf (histo,
"EE E25" ) ;
132 meEEe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
134 sprintf (histo,
"EE E1oE4" ) ;
135 meEEe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
137 sprintf (histo,
"EE E1oE9" ) ;
138 meEEe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
140 sprintf (histo,
"EE E4oE9" ) ;
141 meEEe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
143 sprintf (histo,
"EE E9oE16" ) ;
144 meEEe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
146 sprintf (histo,
"EE E1oE25" ) ;
147 meEEe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
149 sprintf (histo,
"EE E9oE25" ) ;
150 meEEe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
152 sprintf (histo,
"EE E16oE25" ) ;
153 meEEe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
181 if( ! EcalHitsEE.
isValid() )
return;
186 std::vector<PCaloHit> theEECaloHits;
187 theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
191 std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
201 std::vector<double> econtr(140, 0. );
202 std::vector<double> econtr25(140, 0. );
207 uint32_t nEEzpHits = 0;
208 uint32_t nEEzmHits = 0;
210 for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin();
211 isim != theEECaloHits.end(); ++isim){
213 if ( isim->time() > 500. ) {
continue; }
215 CaloHitMap[ isim->id()].push_back(&(*isim));
220 <<
" CaloHit " << isim->getName() <<
"\n"
221 <<
" DetID = "<<isim->id()<<
" EEDetId = " << eeid.ix() <<
" " << eeid.iy() <<
"\n"
222 <<
" Time = " << isim->time() <<
"\n"
223 <<
" Track Id = " << isim->geantTrackId() <<
"\n"
224 <<
" Energy = " << isim->energy();
226 uint32_t crystid = eeid.rawId();
228 if (eeid.zside() > 0 ) {
230 EEetzp_ += isim->energy();
231 eemapzp[crystid] += isim->energy();
234 else if (eeid.zside() < 0 ) {
236 EEetzm_ += isim->energy();
237 eemapzm[crystid] += isim->energy();
242 if( isim->energy() > 0 ) {
244 int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
245 if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
248 eemap[crystid] += isim->energy();
257 for (
std::map<uint32_t,
float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it )
meEEcrystalEnergy_->
Fill((*it).second);
267 int nEEHits = nEEzmHits + nEEzpHits;
272 int bx = myEEid.
ix();
273 int by = myEEid.
iy();
274 int bz = myEEid.
zside();
282 std::vector<uint32_t> ids25; ids25 =
getIdsAroundMax(5,5,bx,by,bz,eemap);
284 for(
unsigned i=0;
i<25;
i++ ) {
285 for(
unsigned int j=0;
j<CaloHitMap[ids25[
i]].size();
j++ ) {
286 if( CaloHitMap[ids25[
i]][
j]->
energy() > 0 ) {
287 int log10i = int( ( log10( CaloHitMap[ids25[
i]][
j]->
energy()) + 10. ) * 10. );
288 if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[
i]][
j]->energy();
312 for(
int i=0;
i<140;
i++ ) {
318 for(
int i=0;
i<140;
i++ ) {
325 if( MyPEcalValidInfo.
isValid() ) {
326 if ( MyPEcalValidInfo->ee1x1() > 0. ) {
327 std::vector<float> EX0 = MyPEcalValidInfo->eX0();
329 for (
int myStep=0; myStep< 26; myStep++ ) {
338 int centralX,
int centralY,
339 int centralZ,
MapType& themap){
342 float totalEnergy = 0.;
344 int goBackInX = nCellInX/2;
345 int goBackInY = nCellInY/2;
346 int startX = centralX-goBackInX;
347 int startY = centralY-goBackInY;
349 for (
int ix=startX; ix<startX+nCellInX; ix++) {
350 for (
int iy=startY; iy<startY+nCellInY; iy++) {
358 totalEnergy += themap[
index];
364 << nCellInX <<
" x " << nCellInY
365 <<
" EE matrix energy = " << totalEnergy
366 <<
" for " << ncristals <<
" crystals";
374 std::vector<uint32_t> ids( nCellInX*nCellInY );
376 int goBackInX = nCellInX/2;
377 int goBackInY = nCellInY/2;
378 int startX = centralX-goBackInX;
379 int startY = centralY-goBackInY;
381 for (
int ix=startX; ix<startX+nCellInX; ix++) {
382 for (
int iy=startY; iy<startY+nCellInY; iy++) {
390 ids[ncristals] =
index;
399 int CentralX,
int CentralY,
int CentralZ,
402 int goBackInX = nCellInX/2;
403 int goBackInY = nCellInY/2;
405 int startX = CentralX - goBackInX;
406 int startY = CentralY - goBackInY;
409 for (
int ix = startX; ix < startX+nCellInX; ix ++ ) {
411 for(
int iy = startY; iy < startY + nCellInY; iy++ ) {
419 fillmap[i++] = themap[
index];
424 if ( fillmap[i/2] == themap[centerid] )
433 float e012 = themap[0]+themap[1]+themap[2];
434 float e036 = themap[0]+themap[3]+themap[6];
435 float e678 = themap[6]+themap[7]+themap[8];
436 float e258 = themap[2]+themap[5]+themap[8];
438 if ( (e012>e678 || e012==e678) && (e036>e258 || e036==e258))
439 return E22=themap[0]+themap[1]+themap[3]+themap[4];
440 else if ( (e012>e678 || e012==e678) && (e036<e258 || e036==e258) )
441 return E22=themap[1]+themap[2]+themap[4]+themap[5];
442 else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
443 return E22=themap[3]+themap[4]+themap[6]+themap[7];
444 else if ( (e012<e678|| e012==e678) && (e036<e258|| e036==e258) )
445 return E22=themap[4]+themap[5]+themap[7]+themap[8];
453 float e0_4 = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
454 float e0_20 = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
455 float e4_24 = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
456 float e0_24 = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
458 if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
459 return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
460 else if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
461 return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
462 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
463 return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
464 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
465 return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
474 uint32_t unitWithMaxEnergy = 0;
475 float maxEnergy = 0.;
477 MapType::iterator iter;
478 for (iter = themap.begin(); iter != themap.end(); iter++) {
480 if (maxEnergy < (*iter).second) {
481 maxEnergy = (*iter).second;
482 unitWithMaxEnergy = (*iter).first;
487 <<
" max energy of " << maxEnergy
488 <<
" GeV in Unit id " << unitWithMaxEnergy;
489 return unitWithMaxEnergy;
EventNumber_t event() const
MonitorElement * meEEe9oe16_
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * meEEzpCrystals_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
edm::EDGetTokenT< PEcalValidInfo > ValidationCollectionToken
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_
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_
edm::EDGetTokenT< edm::PCaloHitContainer > EEHitsToken
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.