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  EBHitsToken = consumes <edm::PCaloHitContainer> (edm::InputTag(std::string(g4InfoLabel),std::string(EBHitsCollection)));
24  // verbosity switch
25  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
26 
27  // get hold of back-end interface
28  dbe_ = 0;
29  dbe_ = edm::Service<DQMStore>().operator->();
30  if ( dbe_ ) {
31  if ( verbose_ ) { dbe_->setVerbose(1); }
32  else { dbe_->setVerbose(0); }
33  }
34 
35  if ( dbe_ ) {
36  if ( verbose_ ) dbe_->showDirStructure();
37  }
38 
39  menEBHits_ = 0;
40  menEBCrystals_ = 0;
41  meEBOccupancy_ = 0;
43  meEBhitEnergy_ = 0;
47  meEBhitEnergy2_ = 0;
50 
51  meEBe1_ = 0;
52  meEBe4_ = 0;
53  meEBe9_ = 0;
54  meEBe16_ = 0;
55  meEBe25_ = 0;
56 
57  meEBe1oe4_ = 0;
58  meEBe1oe9_ = 0;
59  meEBe4oe9_ = 0;
60  meEBe9oe16_ = 0;
61  meEBe1oe25_ = 0;
62  meEBe9oe25_ = 0;
63  meEBe16oe25_ = 0;
64 
65  myEntries = 0;
66  for (int myStep = 0; myStep<26; myStep++) { eRLength[myStep] = 0.0; }
67 
68  Char_t histo[200];
69 
70  if ( dbe_ ) {
71  dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
72 
73  sprintf (histo, "EB hits multiplicity" ) ;
74  menEBHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
75 
76  sprintf (histo, "EB crystals multiplicity" ) ;
77  menEBCrystals_ = dbe_->book1D(histo, histo, 200, 0., 2000.) ;
78 
79  sprintf (histo, "EB occupancy" ) ;
80  meEBOccupancy_ = dbe_->book2D(histo, histo, 360, 0., 360., 170, -85., 85.);
81 
82  sprintf (histo, "EB longitudinal shower profile" ) ;
83  meEBLongitudinalShower_ = dbe_->bookProfile(histo, histo, 26, 0., 26., 100, 0., 20000.);
84 
85  sprintf (histo, "EB hits energy spectrum" );
86  meEBhitEnergy_ = dbe_->book1D(histo, histo, 4000, 0., 400.);
87 
88  sprintf (histo, "EB hits log10energy spectrum" );
89  meEBhitLog10Energy_ = dbe_->book1D(histo, histo, 140, -10., 4.);
90 
91  sprintf (histo, "EB hits log10energy spectrum vs normalized energy" );
92  meEBhitLog10EnergyNorm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
93 
94  sprintf (histo, "EB hits log10energy spectrum vs normalized energy25" );
95  meEBhitLog10Energy25Norm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
96 
97  sprintf (histo, "EB hits energy spectrum 2" );
98  meEBhitEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
99 
100  sprintf (histo, "EB crystal energy spectrum" );
101  meEBcrystalEnergy_ = dbe_->book1D(histo, histo, 5000, 0., 50.);
102 
103  sprintf (histo, "EB crystal energy spectrum 2" );
104  meEBcrystalEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
105 
106  sprintf (histo, "EB E1" ) ;
107  meEBe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
108 
109  sprintf (histo, "EB E4" ) ;
110  meEBe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
111 
112  sprintf (histo, "EB E9" ) ;
113  meEBe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
114 
115  sprintf (histo, "EB E16" ) ;
116  meEBe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
117 
118  sprintf (histo, "EB E25" ) ;
119  meEBe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
120 
121  sprintf (histo, "EB E1oE4" ) ;
122  meEBe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
123 
124  sprintf (histo, "EB E1oE9" ) ;
125  meEBe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
126 
127  sprintf (histo, "EB E4oE9" ) ;
128  meEBe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
129 
130  sprintf (histo, "EB E9oE16" ) ;
131  meEBe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
132 
133  sprintf (histo, "EB E1oE25" ) ;
134  meEBe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
135 
136  sprintf (histo, "EB E9oE25" ) ;
137  meEBe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
138 
139  sprintf (histo, "EB E16oE25" ) ;
140  meEBe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
141  }
142 }
143 
145 
146 }
147 
149 
150 }
151 
153 
154  //for (int myStep=0; myStep<26; myStep++){
155  // if (meEBLongitudinalShower_) meEBLongitudinalShower_->Fill(float(myStep), eRLength[myStep]/myEntries);
156  //}
157 }
158 
160 
161  edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
162 
164  e.getByToken(EBHitsToken,EcalHitsEB);
165 
166  // Do nothing if no Barrel data available
167  if( ! EcalHitsEB.isValid() ) return;
168 
169  edm::Handle<PEcalValidInfo> MyPEcalValidInfo;
170  e.getByToken(ValidationCollectionToken,MyPEcalValidInfo);
171 
172  std::vector<PCaloHit> theEBCaloHits;
173  theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
174 
175  myEntries++;
176 
177  double EBEnergy_ = 0.;
178  std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
179 
180  double eb1 = 0.0;
181  double eb4 = 0.0;
182  double eb9 = 0.0;
183  double eb16 = 0.0;
184  double eb25 = 0.0;
185  std::vector<double> econtr(140, 0. );
186  std::vector<double> econtr25(140, 0. );
187 
188  MapType ebmap;
189  uint32_t nEBHits = 0;
190 
191  for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin();
192  isim != theEBCaloHits.end(); ++isim){
193 
194  if ( isim->time() > 500. ) { continue; }
195 
196  CaloHitMap[ isim->id()].push_back( &(*isim) );
197 
198  EBDetId ebid (isim->id()) ;
199 
200  LogDebug("HitInfo")
201  << " CaloHit " << isim->getName() << "\n"
202  << " DetID = " << isim->id() << " EBDetId = " << ebid.ieta() << " " << ebid.iphi() << "\n"
203  << " Time = " << isim->time() << "\n"
204  << " Track Id = " << isim->geantTrackId() << "\n"
205  << " Energy = " << isim->energy();
206 
207  if (meEBOccupancy_) meEBOccupancy_->Fill( ebid.iphi(), ebid.ieta() );
208 
209  uint32_t crystid = ebid.rawId();
210  ebmap[crystid] += isim->energy();
211 
212  EBEnergy_ += isim->energy();
213  nEBHits++;
214  meEBhitEnergy_->Fill(isim->energy());
215  if( isim->energy() > 0 ) {
216  meEBhitLog10Energy_->Fill(log10(isim->energy()));
217  int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
218  if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
219  }
220  meEBhitEnergy2_->Fill(isim->energy());
221 
222  }
223 
224  if (menEBCrystals_) menEBCrystals_->Fill(ebmap.size());
225  if (meEBcrystalEnergy_) {
226  for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy_->Fill((*it).second);
227  }
228  if (meEBcrystalEnergy2_) {
229  for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy2_->Fill((*it).second);
230  }
231 
232  if (menEBHits_) menEBHits_->Fill(nEBHits);
233 
234  if (nEBHits > 0 ) {
235 
236  uint32_t ebcenterid = getUnitWithMaxEnergy(ebmap);
237  EBDetId myEBid(ebcenterid);
238  int bx = myEBid.ietaAbs();
239  int by = myEBid.iphi();
240  int bz = myEBid.zside();
241  eb1 = energyInMatrixEB(1,1,bx,by,bz,ebmap);
242  if (meEBe1_) meEBe1_->Fill(eb1);
243  eb9 = energyInMatrixEB(3,3,bx,by,bz,ebmap);
244  if (meEBe9_) meEBe9_->Fill(eb9);
245  eb25= energyInMatrixEB(5,5,bx,by,bz,ebmap);
246  if (meEBe25_) meEBe25_->Fill(eb25);
247 
248  std::vector<uint32_t> ids25; ids25 = getIdsAroundMax(5,5,bx,by,bz,ebmap);
249 
250  for( unsigned i=0; i<25; i++ ) {
251  for( unsigned int j=0; j<CaloHitMap[ids25[i]].size(); j++ ) {
252  if( CaloHitMap[ids25[i]][j]->energy() > 0 ) {
253  int log10i = int( ( log10( CaloHitMap[ids25[i]][j]->energy()) + 10. ) * 10. );
254  if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[i]][j]->energy();
255  }
256  }
257  }
258 
259  MapType newebmap;
260  if( fillEBMatrix(3,3,bx,by,bz,newebmap, ebmap)){
261  eb4 = eCluster2x2(newebmap);
262  if (meEBe4_) meEBe4_->Fill(eb4);
263  }
264  if( fillEBMatrix(5,5,bx,by,bz,newebmap, ebmap)){
265  eb16 = eCluster4x4(eb9,newebmap);
266  if (meEBe16_) meEBe16_->Fill(eb16);
267  }
268 
269  if (meEBe1oe4_ && eb4 > 0.1 ) meEBe1oe4_ -> Fill(eb1/eb4);
270  if (meEBe1oe9_ && eb9 > 0.1 ) meEBe1oe9_ -> Fill(eb1/eb9);
271  if (meEBe4oe9_ && eb9 > 0.1 ) meEBe4oe9_ -> Fill(eb4/eb9);
272  if (meEBe9oe16_ && eb16 > 0.1 ) meEBe9oe16_ -> Fill(eb9/eb16);
273  if (meEBe1oe25_ && eb25 > 0.1 ) meEBe1oe25_ -> Fill(eb1/eb25);
274  if (meEBe9oe25_ && eb25 > 0.1 ) meEBe9oe25_ -> Fill(eb9/eb25);
275  if (meEBe16oe25_ && eb25 > 0.1 ) meEBe16oe25_-> Fill(eb16/eb25);
276 
277  if( meEBhitLog10EnergyNorm_ && EBEnergy_ != 0 ) {
278  for( int i=0; i<140; i++ ) {
279  meEBhitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/EBEnergy_ );
280  }
281  }
282 
283  if( meEBhitLog10Energy25Norm_ && eb25 != 0 ) {
284  for( int i=0; i<140; i++ ) {
285  meEBhitLog10Energy25Norm_->Fill( -10.+(float(i)+0.5)/10., econtr25[i]/eb25 );
286  }
287  }
288 
289 
290  }
291 
292 
293  if( MyPEcalValidInfo.isValid() ) {
294  if ( MyPEcalValidInfo->eb1x1() > 0. ) {
295  std::vector<float> BX0 = MyPEcalValidInfo->bX0();
297  for (int myStep=0; myStep< 26; myStep++ ) {
298  eRLength[myStep] += BX0[myStep];
300  }
301  }
302  }
303 }
304 
305 float EcalBarrelSimHitsValidation::energyInMatrixEB(int nCellInEta, int nCellInPhi,
306  int centralEta, int centralPhi,
307  int centralZ, MapType& themap){
308 
309  int ncristals = 0;
310  float totalEnergy = 0.;
311 
312  int goBackInEta = nCellInEta/2;
313  int goBackInPhi = nCellInPhi/2;
314  int startEta = centralZ*centralEta-goBackInEta;
315  int startPhi = centralPhi-goBackInPhi;
316 
317  for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
318  for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
319 
320  uint32_t index ;
321  if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
322  if (iphi< 1) { index = EBDetId(ieta,iphi+360).rawId(); }
323  else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
324  else { index = EBDetId(ieta,iphi).rawId(); }
325 
326  totalEnergy += themap[index];
327  ncristals += 1;
328  }
329  }
330 
331  LogDebug("GeomInfo")
332  << nCellInEta << " x " << nCellInPhi
333  << " EB matrix energy = " << totalEnergy
334  << " for " << ncristals << " crystals" ;
335  return totalEnergy;
336 
337 }
338 
339 std::vector<uint32_t> EcalBarrelSimHitsValidation::getIdsAroundMax( int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType& themap){
340 
341  int ncristals = 0;
342  std::vector<uint32_t> ids( nCellInEta*nCellInPhi );
343 
344  int goBackInEta = nCellInEta/2;
345  int goBackInPhi = nCellInPhi/2;
346  int startEta = centralZ*centralEta-goBackInEta;
347  int startPhi = centralPhi-goBackInPhi;
348 
349  for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
350  for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
351 
352  uint32_t index ;
353  if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
354  if (iphi< 1) { index = EBDetId(ieta,iphi+360).rawId(); }
355  else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
356  else { index = EBDetId(ieta,iphi).rawId(); }
357  ids[ncristals] = index;
358  ncristals += 1;
359  }
360  }
361 
362  return ids;
363 
364 }
365 
366 bool EcalBarrelSimHitsValidation::fillEBMatrix(int nCellInEta, int nCellInPhi,
367  int CentralEta, int CentralPhi,int CentralZ,
368  MapType& fillmap, MapType& themap)
369 {
370  int goBackInEta = nCellInEta/2;
371  int goBackInPhi = nCellInPhi/2;
372 
373  int startEta = CentralZ*CentralEta - goBackInEta;
374  int startPhi = CentralPhi - goBackInPhi;
375 
376  int i = 0 ;
377  for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
378  for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
379  uint32_t index;
380  if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
381  if (iphi< 1) { index = EBDetId(ieta,iphi+360).rawId(); }
382  else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
383  else { index = EBDetId(ieta,iphi).rawId(); }
384  fillmap[i++] = themap[index];
385  }
386  }
387 
388  uint32_t ebcenterid = getUnitWithMaxEnergy(themap);
389 
390  if ( fillmap[i/2] == themap[ebcenterid] )
391  return true;
392  else
393  return false;
394 }
395 
397  float E22=0.;
398  float e012 = themap[0]+themap[1]+themap[2];
399  float e036 = themap[0]+themap[3]+themap[6];
400  float e678 = themap[6]+themap[7]+themap[8];
401  float e258 = themap[2]+themap[5]+themap[8];
402 
403  if ( (e012>e678 || e012==e678) && (e036>e258 || e036==e258))
404  return E22=themap[0]+themap[1]+themap[3]+themap[4];
405  else if ( (e012>e678 || e012==e678) && (e036<e258 || e036==e258) )
406  return E22=themap[1]+themap[2]+themap[4]+themap[5];
407  else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
408  return E22=themap[3]+themap[4]+themap[6]+themap[7];
409  else if ( (e012<e678|| e012==e678) && (e036<e258|| e036==e258) )
410  return E22=themap[4]+themap[5]+themap[7]+themap[8];
411  else {
412  return E22;
413  }
414 }
415 
417  float E44=0.;
418  float e0_4 = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
419  float e0_20 = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
420  float e4_24 = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
421  float e0_24 = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
422 
423  if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
424  return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
425  else if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
426  return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
427  else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
428  return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
429  else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
430  return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
431  else{
432  return E44;
433  }
434 }
435 
437 
438  //look for max
439  uint32_t unitWithMaxEnergy = 0;
440  float maxEnergy = 0.;
441 
442  MapType::iterator iter;
443  for (iter = themap.begin(); iter != themap.end(); iter++) {
444 
445  if (maxEnergy < (*iter).second) {
446  maxEnergy = (*iter).second;
447  unitWithMaxEnergy = (*iter).first;
448  }
449  }
450 
451  LogDebug("GeomInfo")
452  << " max energy of " << maxEnergy
453  << " GeV in Unit id " << unitWithMaxEnergy;
454  return unitWithMaxEnergy;
455 }
456 
457 
#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)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
void analyze(const edm::Event &e, const edm::EventSetup &c)
Analyze.
void Fill(long long x)
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
float eCluster4x4(float e33, MapType &themap)
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
edm::EDGetTokenT< edm::PCaloHitContainer > EBHitsToken
bool isValid() const
Definition: HandleBase.h:76
edm::EDGetTokenT< PEcalValidInfo > ValidationCollectionToken
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:49
int zside() const
get the z-side of the crystal (1/-1)
Definition: EBDetId.h:47
EcalBarrelSimHitsValidation(const edm::ParameterSet &ps)
Constructor.