CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalBarrelSimHitsValidation.cc
Go to the documentation of this file.
1 /*
2  * \file EcalBarrelSimHitsValidation.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  EBHitsCollection(ps.getParameter<std::string>("EBHitsCollection")),
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  menEBHits_ = 0;
38  menEBCrystals_ = 0;
39  meEBOccupancy_ = 0;
41  meEBhitEnergy_ = 0;
45  meEBhitEnergy2_ = 0;
48 
49  meEBe1_ = 0;
50  meEBe4_ = 0;
51  meEBe9_ = 0;
52  meEBe16_ = 0;
53  meEBe25_ = 0;
54 
55  meEBe1oe4_ = 0;
56  meEBe1oe9_ = 0;
57  meEBe4oe9_ = 0;
58  meEBe9oe16_ = 0;
59  meEBe1oe25_ = 0;
60  meEBe9oe25_ = 0;
61  meEBe16oe25_ = 0;
62 
63  myEntries = 0;
64  for (int myStep = 0; myStep<26; myStep++) { eRLength[myStep] = 0.0; }
65 
66  Char_t histo[200];
67 
68  if ( dbe_ ) {
69  dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
70 
71  sprintf (histo, "EB hits multiplicity" ) ;
72  menEBHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
73 
74  sprintf (histo, "EB crystals multiplicity" ) ;
75  menEBCrystals_ = dbe_->book1D(histo, histo, 200, 0., 2000.) ;
76 
77  sprintf (histo, "EB occupancy" ) ;
78  meEBOccupancy_ = dbe_->book2D(histo, histo, 360, 0., 360., 170, -85., 85.);
79 
80  sprintf (histo, "EB longitudinal shower profile" ) ;
81  meEBLongitudinalShower_ = dbe_->bookProfile(histo, histo, 26, 0., 26., 100, 0., 20000.);
82 
83  sprintf (histo, "EB hits energy spectrum" );
84  meEBhitEnergy_ = dbe_->book1D(histo, histo, 4000, 0., 400.);
85 
86  sprintf (histo, "EB hits log10energy spectrum" );
87  meEBhitLog10Energy_ = dbe_->book1D(histo, histo, 140, -10., 4.);
88 
89  sprintf (histo, "EB hits log10energy spectrum vs normalized energy" );
90  meEBhitLog10EnergyNorm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
91 
92  sprintf (histo, "EB hits log10energy spectrum vs normalized energy25" );
93  meEBhitLog10Energy25Norm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
94 
95  sprintf (histo, "EB hits energy spectrum 2" );
96  meEBhitEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
97 
98  sprintf (histo, "EB crystal energy spectrum" );
99  meEBcrystalEnergy_ = dbe_->book1D(histo, histo, 5000, 0., 50.);
100 
101  sprintf (histo, "EB crystal energy spectrum 2" );
102  meEBcrystalEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
103 
104  sprintf (histo, "EB E1" ) ;
105  meEBe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
106 
107  sprintf (histo, "EB E4" ) ;
108  meEBe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
109 
110  sprintf (histo, "EB E9" ) ;
111  meEBe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
112 
113  sprintf (histo, "EB E16" ) ;
114  meEBe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
115 
116  sprintf (histo, "EB E25" ) ;
117  meEBe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
118 
119  sprintf (histo, "EB E1oE4" ) ;
120  meEBe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
121 
122  sprintf (histo, "EB E1oE9" ) ;
123  meEBe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
124 
125  sprintf (histo, "EB E4oE9" ) ;
126  meEBe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
127 
128  sprintf (histo, "EB E9oE16" ) ;
129  meEBe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
130 
131  sprintf (histo, "EB E1oE25" ) ;
132  meEBe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
133 
134  sprintf (histo, "EB E9oE25" ) ;
135  meEBe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
136 
137  sprintf (histo, "EB E16oE25" ) ;
138  meEBe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
139  }
140 }
141 
143 
144 }
145 
147 
148 }
149 
151 
152  //for (int myStep=0; myStep<26; myStep++){
153  // if (meEBLongitudinalShower_) meEBLongitudinalShower_->Fill(float(myStep), eRLength[myStep]/myEntries);
154  //}
155 }
156 
158 
159  edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
160 
163 
164  // Do nothing if no Barrel data available
165  if( ! EcalHitsEB.isValid() ) return;
166 
167  edm::Handle<PEcalValidInfo> MyPEcalValidInfo;
168  e.getByLabel(g4InfoLabel,ValidationCollection,MyPEcalValidInfo);
169 
170  std::vector<PCaloHit> theEBCaloHits;
171  theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
172 
173  myEntries++;
174 
175  double EBEnergy_ = 0.;
176  std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
177 
178  double eb1 = 0.0;
179  double eb4 = 0.0;
180  double eb9 = 0.0;
181  double eb16 = 0.0;
182  double eb25 = 0.0;
183  std::vector<double> econtr(140, 0. );
184  std::vector<double> econtr25(140, 0. );
185 
186  MapType ebmap;
187  uint32_t nEBHits = 0;
188 
189  for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin();
190  isim != theEBCaloHits.end(); ++isim){
191 
192  if ( isim->time() > 500. ) { continue; }
193 
194  CaloHitMap[ isim->id()].push_back( &(*isim) );
195 
196  EBDetId ebid (isim->id()) ;
197 
198  LogDebug("HitInfo")
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();
204 
205  if (meEBOccupancy_) meEBOccupancy_->Fill( ebid.iphi(), ebid.ieta() );
206 
207  uint32_t crystid = ebid.rawId();
208  ebmap[crystid] += isim->energy();
209 
210  EBEnergy_ += isim->energy();
211  nEBHits++;
212  meEBhitEnergy_->Fill(isim->energy());
213  if( isim->energy() > 0 ) {
214  meEBhitLog10Energy_->Fill(log10(isim->energy()));
215  int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
216  if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
217  }
218  meEBhitEnergy2_->Fill(isim->energy());
219 
220  }
221 
222  if (menEBCrystals_) menEBCrystals_->Fill(ebmap.size());
223  if (meEBcrystalEnergy_) {
224  for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy_->Fill((*it).second);
225  }
226  if (meEBcrystalEnergy2_) {
227  for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy2_->Fill((*it).second);
228  }
229 
230  if (menEBHits_) menEBHits_->Fill(nEBHits);
231 
232  if (nEBHits > 0 ) {
233 
234  uint32_t ebcenterid = getUnitWithMaxEnergy(ebmap);
235  EBDetId myEBid(ebcenterid);
236  int bx = myEBid.ietaAbs();
237  int by = myEBid.iphi();
238  int bz = myEBid.zside();
239  eb1 = energyInMatrixEB(1,1,bx,by,bz,ebmap);
240  if (meEBe1_) meEBe1_->Fill(eb1);
241  eb9 = energyInMatrixEB(3,3,bx,by,bz,ebmap);
242  if (meEBe9_) meEBe9_->Fill(eb9);
243  eb25= energyInMatrixEB(5,5,bx,by,bz,ebmap);
244  if (meEBe25_) meEBe25_->Fill(eb25);
245 
246  std::vector<uint32_t> ids25; ids25 = getIdsAroundMax(5,5,bx,by,bz,ebmap);
247 
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();
253  }
254  }
255  }
256 
257  MapType newebmap;
258  if( fillEBMatrix(3,3,bx,by,bz,newebmap, ebmap)){
259  eb4 = eCluster2x2(newebmap);
260  if (meEBe4_) meEBe4_->Fill(eb4);
261  }
262  if( fillEBMatrix(5,5,bx,by,bz,newebmap, ebmap)){
263  eb16 = eCluster4x4(eb9,newebmap);
264  if (meEBe16_) meEBe16_->Fill(eb16);
265  }
266 
267  if (meEBe1oe4_ && eb4 > 0.1 ) meEBe1oe4_ -> Fill(eb1/eb4);
268  if (meEBe1oe9_ && eb9 > 0.1 ) meEBe1oe9_ -> Fill(eb1/eb9);
269  if (meEBe4oe9_ && eb9 > 0.1 ) meEBe4oe9_ -> Fill(eb4/eb9);
270  if (meEBe9oe16_ && eb16 > 0.1 ) meEBe9oe16_ -> Fill(eb9/eb16);
271  if (meEBe1oe25_ && eb25 > 0.1 ) meEBe1oe25_ -> Fill(eb1/eb25);
272  if (meEBe9oe25_ && eb25 > 0.1 ) meEBe9oe25_ -> Fill(eb9/eb25);
273  if (meEBe16oe25_ && eb25 > 0.1 ) meEBe16oe25_-> Fill(eb16/eb25);
274 
275  if( meEBhitLog10EnergyNorm_ && EBEnergy_ != 0 ) {
276  for( int i=0; i<140; i++ ) {
277  meEBhitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/EBEnergy_ );
278  }
279  }
280 
281  if( meEBhitLog10Energy25Norm_ && eb25 != 0 ) {
282  for( int i=0; i<140; i++ ) {
283  meEBhitLog10Energy25Norm_->Fill( -10.+(float(i)+0.5)/10., econtr25[i]/eb25 );
284  }
285  }
286 
287 
288  }
289 
290 
291  if( MyPEcalValidInfo.isValid() ) {
292  if ( MyPEcalValidInfo->eb1x1() > 0. ) {
293  std::vector<float> BX0 = MyPEcalValidInfo->bX0();
295  for (int myStep=0; myStep< 26; myStep++ ) {
296  eRLength[myStep] += BX0[myStep];
298  }
299  }
300  }
301 }
302 
303 float EcalBarrelSimHitsValidation::energyInMatrixEB(int nCellInEta, int nCellInPhi,
304  int centralEta, int centralPhi,
305  int centralZ, MapType& themap){
306 
307  int ncristals = 0;
308  float totalEnergy = 0.;
309 
310  int goBackInEta = nCellInEta/2;
311  int goBackInPhi = nCellInPhi/2;
312  int startEta = centralZ*centralEta-goBackInEta;
313  int startPhi = centralPhi-goBackInPhi;
314 
315  for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
316  for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
317 
318  uint32_t index ;
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(); }
322  else { index = EBDetId(ieta,iphi).rawId(); }
323 
324  totalEnergy += themap[index];
325  ncristals += 1;
326  }
327  }
328 
329  LogDebug("GeomInfo")
330  << nCellInEta << " x " << nCellInPhi
331  << " EB matrix energy = " << totalEnergy
332  << " for " << ncristals << " crystals" ;
333  return totalEnergy;
334 
335 }
336 
337 std::vector<uint32_t> EcalBarrelSimHitsValidation::getIdsAroundMax( int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType& themap){
338 
339  int ncristals = 0;
340  std::vector<uint32_t> ids( nCellInEta*nCellInPhi );
341 
342  int goBackInEta = nCellInEta/2;
343  int goBackInPhi = nCellInPhi/2;
344  int startEta = centralZ*centralEta-goBackInEta;
345  int startPhi = centralPhi-goBackInPhi;
346 
347  for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
348  for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
349 
350  uint32_t index ;
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(); }
354  else { index = EBDetId(ieta,iphi).rawId(); }
355  ids[ncristals] = index;
356  ncristals += 1;
357  }
358  }
359 
360  return ids;
361 
362 }
363 
364 bool EcalBarrelSimHitsValidation::fillEBMatrix(int nCellInEta, int nCellInPhi,
365  int CentralEta, int CentralPhi,int CentralZ,
366  MapType& fillmap, MapType& themap)
367 {
368  int goBackInEta = nCellInEta/2;
369  int goBackInPhi = nCellInPhi/2;
370 
371  int startEta = CentralZ*CentralEta - goBackInEta;
372  int startPhi = CentralPhi - goBackInPhi;
373 
374  int i = 0 ;
375  for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
376  for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
377  uint32_t index;
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(); }
381  else { index = EBDetId(ieta,iphi).rawId(); }
382  fillmap[i++] = themap[index];
383  }
384  }
385 
386  uint32_t ebcenterid = getUnitWithMaxEnergy(themap);
387 
388  if ( fillmap[i/2] == themap[ebcenterid] )
389  return true;
390  else
391  return false;
392 }
393 
395  float E22=0.;
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];
400 
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];
409  else {
410  return E22;
411  }
412 }
413 
415  float E44=0.;
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];
420 
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];
429  else{
430  return E44;
431  }
432 }
433 
435 
436  //look for max
437  uint32_t unitWithMaxEnergy = 0;
438  float maxEnergy = 0.;
439 
440  MapType::iterator iter;
441  for (iter = themap.begin(); iter != themap.end(); iter++) {
442 
443  if (maxEnergy < (*iter).second) {
444  maxEnergy = (*iter).second;
445  unitWithMaxEnergy = (*iter).first;
446  }
447  }
448 
449  LogDebug("GeomInfo")
450  << " max energy of " << maxEnergy
451  << " GeV in Unit id " << unitWithMaxEnergy;
452  return unitWithMaxEnergy;
453 }
454 
455 
#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
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)
void analyze(const edm::Event &e, const edm::EventSetup &c)
Analyze.
#define abs(x)
Definition: mlp_lapack.h:159
void Fill(long long x)
int iphi() const
get the crystal iphi
Definition: EBDetId.h:46
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
float eCluster4x4(float e33, MapType &themap)
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
int j
Definition: DBlmapReader.cc:9
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
uint32_t getUnitWithMaxEnergy(MapType &themap)
virtual float energyInMatrixEB(int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType &themap)
edm::EventID id() const
Definition: EventBase.h:56
bool fillEBMatrix(int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &fillmap, MapType &themap)
void Reset(void)
reset ME (ie. contents, errors, etc)
int ietaAbs() const
get the absolute value of the crystal ieta
Definition: EBDetId.h:42
int zside() const
get the z-side of the crystal (1/-1)
Definition: EBDetId.h:40
EcalBarrelSimHitsValidation(const edm::ParameterSet &ps)
Constructor.