CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalEndcapSimHitsValidation.cc
Go to the documentation of this file.
1 /*
2  * \file EcalEndcapSimHitsValidation.cc
3  *
4  * \author C.Rovelli
5  *
6 */
7 
12 
13 using namespace cms;
14 using namespace edm;
15 using namespace std;
16 
18  g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
19  EEHitsCollection(ps.getParameter<std::string>("EEHitsCollection")),
20  ValidationCollection(ps.getParameter<std::string>("ValidationCollection")){
21 
22  // verbosity switch
23  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
24 
25  // get hold of back-end interface
26  dbe_ = 0;
27  dbe_ = edm::Service<DQMStore>().operator->();
28  if ( dbe_ ) {
29  if ( verbose_ ) { dbe_->setVerbose(1); }
30  else { dbe_->setVerbose(0); }
31  }
32 
33  if ( dbe_ ) {
34  if ( verbose_ ) dbe_->showDirStructure();
35  }
36 
37 
38  meEEzpHits_ = 0;
39  meEEzmHits_ = 0;
40  meEEzpCrystals_ = 0;
41  meEEzmCrystals_ = 0;
42  meEEzpOccupancy_ = 0;
43  meEEzmOccupancy_ = 0;
45  meEEHitEnergy_ = 0;
49  meEEHitEnergy2_ = 0;
52 
53  meEEe1_ = 0;
54  meEEe4_ = 0;
55  meEEe9_ = 0;
56  meEEe16_ = 0;
57  meEEe25_ = 0;
58 
59  meEEe1oe4_ = 0;
60  meEEe1oe9_ = 0;
61  meEEe4oe9_ = 0;
62  meEEe9oe16_ = 0;
63  meEEe1oe25_ = 0;
64  meEEe9oe25_ = 0;
65  meEEe16oe25_ = 0;
66 
67  myEntries = 0;
68  for ( int myStep = 0; myStep<26; myStep++) { eRLength[myStep] = 0.0; }
69 
70  Char_t histo[200];
71 
72  if ( dbe_ ) {
73  dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
74 
75  sprintf (histo, "EE+ hits multiplicity" ) ;
76  meEEzpHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
77 
78  sprintf (histo, "EE- hits multiplicity" ) ;
79  meEEzmHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
80 
81  sprintf (histo, "EE+ crystals multiplicity" ) ;
82  meEEzpCrystals_ = dbe_->book1D(histo, histo, 200, 0., 2000.) ;
83 
84  sprintf (histo, "EE- crystals multiplicity" ) ;
85  meEEzmCrystals_ = dbe_->book1D(histo, histo, 200, 0., 2000.) ;
86 
87  sprintf (histo, "EE+ occupancy" ) ;
88  meEEzpOccupancy_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
89 
90  sprintf (histo, "EE- occupancy" ) ;
91  meEEzmOccupancy_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
92 
93  sprintf (histo, "EE longitudinal shower profile" ) ;
94  meEELongitudinalShower_ = dbe_->bookProfile(histo, histo, 26,0,26, 100, 0, 20000);
95 
96  sprintf (histo, "EE hits energy spectrum" );
97  meEEHitEnergy_ = dbe_->book1D(histo, histo, 4000, 0., 400.);
98 
99  sprintf (histo, "EE hits log10energy spectrum" );
100  meEEhitLog10Energy_ = dbe_->book1D(histo, histo, 140, -10., 4.);
101 
102  sprintf (histo, "EE hits log10energy spectrum vs normalized energy" );
103  meEEhitLog10EnergyNorm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
104 
105  sprintf (histo, "EE hits log10energy spectrum vs normalized energy25" );
106  meEEhitLog10Energy25Norm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
107 
108  sprintf (histo, "EE hits energy spectrum 2" );
109  meEEHitEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
110 
111  sprintf (histo, "EE crystal energy spectrum" );
112  meEEcrystalEnergy_ = dbe_->book1D(histo, histo, 5000, 0., 50.);
113 
114  sprintf (histo, "EE crystal energy spectrum 2" );
115  meEEcrystalEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
116 
117  sprintf (histo, "EE E1" ) ;
118  meEEe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
119 
120  sprintf (histo, "EE E4" ) ;
121  meEEe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
122 
123  sprintf (histo, "EE E9" ) ;
124  meEEe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
125 
126  sprintf (histo, "EE E16" ) ;
127  meEEe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
128 
129  sprintf (histo, "EE E25" ) ;
130  meEEe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
131 
132  sprintf (histo, "EE E1oE4" ) ;
133  meEEe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
134 
135  sprintf (histo, "EE E1oE9" ) ;
136  meEEe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
137 
138  sprintf (histo, "EE E4oE9" ) ;
139  meEEe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
140 
141  sprintf (histo, "EE E9oE16" ) ;
142  meEEe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
143 
144  sprintf (histo, "EE E1oE25" ) ;
145  meEEe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
146 
147  sprintf (histo, "EE E9oE25" ) ;
148  meEEe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
149 
150  sprintf (histo, "EE E16oE25" ) ;
151  meEEe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
152  }
153 }
154 
156 
157 }
158 
160 
161 }
162 
164 
165  //for ( int myStep = 0; myStep<26; myStep++){
166  // if (meEELongitudinalShower_) meEELongitudinalShower_->Fill(float(myStep), eRLength[myStep]/myEntries);
167  //}
168 
169 }
170 
172 
173  edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
174 
177 
178  // Do nothing if no EndCap data available
179  if( ! EcalHitsEE.isValid() ) return;
180 
181  edm::Handle<PEcalValidInfo> MyPEcalValidInfo;
182  e.getByLabel(g4InfoLabel,ValidationCollection,MyPEcalValidInfo);
183 
184  std::vector<PCaloHit> theEECaloHits;
185  theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
186 
187  myEntries++;
188 
189  std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
190 
191  double EEetzp_ = 0.;
192  double EEetzm_ = 0.;
193 
194  double ee1 = 0.0;
195  double ee4 = 0.0;
196  double ee9 = 0.0;
197  double ee16 = 0.0;
198  double ee25 = 0.0;
199  std::vector<double> econtr(140, 0. );
200  std::vector<double> econtr25(140, 0. );
201 
202  MapType eemap;
203  MapType eemapzp;
204  MapType eemapzm;
205  uint32_t nEEzpHits = 0;
206  uint32_t nEEzmHits = 0;
207 
208  for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin();
209  isim != theEECaloHits.end(); ++isim){
210 
211  if ( isim->time() > 500. ) { continue; }
212 
213  CaloHitMap[ isim->id()].push_back(&(*isim));
214 
215  EEDetId eeid (isim->id()) ;
216 
217  LogDebug("HitInfo")
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();
223 
224  uint32_t crystid = eeid.rawId();
225 
226  if (eeid.zside() > 0 ) {
227  nEEzpHits++;
228  EEetzp_ += isim->energy();
229  eemapzp[crystid] += isim->energy();
230  if (meEEzpOccupancy_) meEEzpOccupancy_->Fill(eeid.ix(), eeid.iy());
231  }
232  else if (eeid.zside() < 0 ) {
233  nEEzmHits++;
234  EEetzm_ += isim->energy();
235  eemapzm[crystid] += isim->energy();
236  if (meEEzmOccupancy_) meEEzmOccupancy_->Fill(eeid.ix(), eeid.iy());
237  }
238 
239  if (meEEHitEnergy_) meEEHitEnergy_->Fill(isim->energy());
240  if( isim->energy() > 0 ) {
241  if( meEEhitLog10Energy_ ) meEEhitLog10Energy_->Fill(log10(isim->energy()));
242  int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
243  if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
244  }
245  if (meEEHitEnergy2_) meEEHitEnergy2_->Fill(isim->energy());
246  eemap[crystid] += isim->energy();
247 
248 
249  }
250 
251  if (meEEzpCrystals_) meEEzpCrystals_->Fill(eemapzp.size());
252  if (meEEzmCrystals_) meEEzmCrystals_->Fill(eemapzm.size());
253 
254  if (meEEcrystalEnergy_) {
255  for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it ) meEEcrystalEnergy_->Fill((*it).second);
256  }
257  if (meEEcrystalEnergy2_) {
258  for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it ) meEEcrystalEnergy2_->Fill((*it).second);
259  }
260 
261  if (meEEzpHits_) meEEzpHits_->Fill(nEEzpHits);
262  if (meEEzmHits_) meEEzmHits_->Fill(nEEzmHits);
263 
264 
265  int nEEHits = nEEzmHits + nEEzpHits;
266  if (nEEHits > 0) {
267 
268  uint32_t eecenterid = getUnitWithMaxEnergy(eemap);
269  EEDetId myEEid(eecenterid);
270  int bx = myEEid.ix();
271  int by = myEEid.iy();
272  int bz = myEEid.zside();
273  ee1 = energyInMatrixEE(1,1,bx,by,bz,eemap);
274  if (meEEe1_) meEEe1_->Fill(ee1);
275  ee9 = energyInMatrixEE(3,3,bx,by,bz,eemap);
276  if (meEEe9_) meEEe9_->Fill(ee9);
277  ee25= energyInMatrixEE(5,5,bx,by,bz,eemap);
278  if (meEEe25_) meEEe25_->Fill(ee25);
279 
280  std::vector<uint32_t> ids25; ids25 = getIdsAroundMax(5,5,bx,by,bz,eemap);
281 
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();
287  }
288  }
289  }
290 
291  MapType neweemap;
292  if( fillEEMatrix(3,3,bx,by,bz,neweemap, eemap)){
293  ee4 = eCluster2x2(neweemap);
294  if (meEEe4_) meEEe4_->Fill(ee4);
295  }
296  if( fillEEMatrix(5,5,bx,by,bz,neweemap, eemap)){
297  ee16 = eCluster4x4(ee9,neweemap);
298  if (meEEe16_) meEEe16_->Fill(ee16);
299  }
300 
301  if (meEEe1oe4_ && ee4 > 0.1 ) meEEe1oe4_ ->Fill(ee1/ee4);
302  if (meEEe1oe9_ && ee9 > 0.1 ) meEEe1oe9_ ->Fill(ee1/ee9);
303  if (meEEe4oe9_ && ee9 > 0.1 ) meEEe4oe9_ ->Fill(ee4/ee9);
304  if (meEEe9oe16_ && ee16 > 0.1 ) meEEe9oe16_ ->Fill(ee9/ee16);
305  if (meEEe16oe25_ && ee25 > 0.1 ) meEEe16oe25_->Fill(ee16/ee25);
306  if (meEEe1oe25_ && ee25 > 0.1 ) meEEe1oe25_ ->Fill(ee1/ee25);
307  if (meEEe9oe25_ && ee25 > 0.1 ) meEEe9oe25_ ->Fill(ee9/ee25);
308 
309  if( meEEhitLog10EnergyNorm_ && (EEetzp_+EEetzm_) != 0 ) {
310  for( int i=0; i<140; i++ ) {
311  meEEhitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/(EEetzp_+EEetzm_) );
312  }
313  }
314 
315  if( meEEhitLog10Energy25Norm_ && ee25 != 0 ) {
316  for( int i=0; i<140; i++ ) {
317  meEEhitLog10Energy25Norm_->Fill( -10.+(float(i)+0.5)/10., econtr25[i]/ee25 );
318  }
319  }
320 
321  }
322 
323  if( MyPEcalValidInfo.isValid() ) {
324  if ( MyPEcalValidInfo->ee1x1() > 0. ) {
325  std::vector<float> EX0 = MyPEcalValidInfo->eX0();
327  for (int myStep=0; myStep< 26; myStep++ ) {
328  eRLength[myStep] += EX0[myStep];
330  }
331  }
332  }
333 }
334 
335 float EcalEndcapSimHitsValidation::energyInMatrixEE(int nCellInX, int nCellInY,
336  int centralX, int centralY,
337  int centralZ, MapType& themap){
338 
339  int ncristals = 0;
340  float totalEnergy = 0.;
341 
342  int goBackInX = nCellInX/2;
343  int goBackInY = nCellInY/2;
344  int startX = centralX-goBackInX;
345  int startY = centralY-goBackInY;
346 
347  for (int ix=startX; ix<startX+nCellInX; ix++) {
348  for (int iy=startY; iy<startY+nCellInY; iy++) {
349  uint32_t index ;
350 
351  if(EEDetId::validDetId(ix,iy,centralZ)) {
352  index = EEDetId(ix,iy,centralZ).rawId();
353  }
354  else { continue; }
355 
356  totalEnergy += themap[index];
357  ncristals += 1;
358  }
359  }
360 
361  LogDebug("GeomInfo")
362  << nCellInX << " x " << nCellInY
363  << " EE matrix energy = " << totalEnergy
364  << " for " << ncristals << " crystals";
365  return totalEnergy;
366 
367 }
368 
369 std::vector<uint32_t> EcalEndcapSimHitsValidation::getIdsAroundMax( int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType& themap){
370 
371  int ncristals = 0;
372  std::vector<uint32_t> ids( nCellInX*nCellInY );
373 
374  int goBackInX = nCellInX/2;
375  int goBackInY = nCellInY/2;
376  int startX = centralX-goBackInX;
377  int startY = centralY-goBackInY;
378 
379  for (int ix=startX; ix<startX+nCellInX; ix++) {
380  for (int iy=startY; iy<startY+nCellInY; iy++) {
381  uint32_t index ;
382 
383  if(EEDetId::validDetId(ix,iy,centralZ)) {
384  index = EEDetId(ix,iy,centralZ).rawId();
385  }
386  else { continue; }
387 
388  ids[ncristals] = index;
389  ncristals += 1;
390  }
391  }
392 
393  return ids;
394 }
395 
396 bool EcalEndcapSimHitsValidation::fillEEMatrix(int nCellInX, int nCellInY,
397  int CentralX, int CentralY,int CentralZ,
398  MapType& fillmap, MapType& themap)
399 {
400  int goBackInX = nCellInX/2;
401  int goBackInY = nCellInY/2;
402 
403  int startX = CentralX - goBackInX;
404  int startY = CentralY - goBackInY;
405 
406  int i = 0 ;
407  for ( int ix = startX; ix < startX+nCellInX; ix ++ ) {
408 
409  for( int iy = startY; iy < startY + nCellInY; iy++ ) {
410 
411  uint32_t index ;
412 
413  if(EEDetId::validDetId(ix,iy,CentralZ)) {
414  index = EEDetId(ix,iy,CentralZ).rawId();
415  }
416  else { continue; }
417  fillmap[i++] = themap[index];
418  }
419  }
420  uint32_t centerid = getUnitWithMaxEnergy(themap);
421 
422  if ( fillmap[i/2] == themap[centerid] )
423  return true;
424  else
425  return false;
426 }
427 
428 
430  float E22=0.;
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];
435 
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];
444  else {
445  return E22;
446  }
447 }
448 
450  float E44=0.;
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];
455 
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];
464  else{
465  return E44;
466  }
467 }
468 
470 
471  //look for max
472  uint32_t unitWithMaxEnergy = 0;
473  float maxEnergy = 0.;
474 
475  MapType::iterator iter;
476  for (iter = themap.begin(); iter != themap.end(); iter++) {
477 
478  if (maxEnergy < (*iter).second) {
479  maxEnergy = (*iter).second;
480  unitWithMaxEnergy = (*iter).first;
481  }
482  }
483 
484  LogDebug("GeomInfo")
485  << " max energy of " << maxEnergy
486  << " GeV in Unit id " << unitWithMaxEnergy;
487  return unitWithMaxEnergy;
488 }
489 
490 
#define LogDebug(id)
RunNumber_t run() const
Definition: EventID.h:42
EventNumber_t event() const
Definition: EventID.h:44
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
int ix() const
Definition: EEDetId.h:76
virtual float energyInMatrixEE(int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap)
void Fill(long long x)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
float eCluster4x4(float e33, MapType &themap)
int zside() const
Definition: EEDetId.h:70
int j
Definition: DBlmapReader.cc:9
int iy() const
Definition: EEDetId.h:82
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:390
bool fillEEMatrix(int nCellInX, int nCellInY, int CentralX, int CentralY, int CentralZ, MapType &fillmap, MapType &themap)
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:248
edm::EventID id() const
Definition: EventBase.h:56
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
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)
EcalEndcapSimHitsValidation(const edm::ParameterSet &ps)
Constructor.