CMS 3D CMS Logo

ECALRecHitAnalyzer.cc
Go to the documentation of this file.
4 
5 // author: Bobby Scurlock, University of Florida
6 // first version 11/20/2006
7 
8 #define DEBUG(X) { if (debug_) { std::cout << X << std::endl; } }
9 
11 {
12 
13  // Retrieve Information from the Configuration File
14  EBRecHitsLabel_ = consumes<EBRecHitCollection> (iConfig.getParameter<edm::InputTag>("EBRecHitsLabel"));
15  EERecHitsLabel_ = consumes<EERecHitCollection> (iConfig.getParameter<edm::InputTag>("EERecHitsLabel"));
16  FolderName_ = iConfig.getUntrackedParameter<std::string>("FolderName");
17  debug_ = iConfig.getParameter<bool>("Debug");
18  // EBRecHitsLabel_= consumes<EcalRecHitCollection>(edm::InputTag(EBRecHitsLabel_));
19  // EERecHitsLabel_= consumes<EcalRecHitCollection>(edm::InputTag(EERecHitsLabel_));
20 }
21 
23  CurrentEvent = -1;
24  // Fill the geometry histograms
25  FillGeometry(iSetup);
26 }
27 
29  edm::Run const & iRun,
30  edm::EventSetup const & )
31 {
32  // get ahold of back-end interface
33  // Book Geometry Histograms
34  ibooker.setCurrentFolder(FolderName_+"/geometry");
35 
36 
37  // ECAL barrel
38  hEB_ieta_iphi_etaMap = ibooker.book2D("hEB_ieta_iphi_etaMap","", 171, -85, 86, 360, 1, 361);
39  hEB_ieta_iphi_phiMap = ibooker.book2D("hEB_ieta_iphi_phiMap","", 171, -85, 86, 360, 1, 361);
40  hEB_ieta_detaMap = ibooker.book1D("hEB_ieta_detaMap","", 171, -85, 86);
41  hEB_ieta_dphiMap = ibooker.book1D("hEB_ieta_dphiMap","", 171, -85, 86);
42  // ECAL +endcap
43  hEEpZ_ix_iy_irMap = ibooker.book2D("hEEpZ_ix_iy_irMap","", 100,1,101, 100,1,101);
44  hEEpZ_ix_iy_xMap = ibooker.book2D("hEEpZ_ix_iy_xMap","", 100,1,101, 100,1,101);
45  hEEpZ_ix_iy_yMap = ibooker.book2D("hEEpZ_ix_iy_yMap","", 100,1,101, 100,1,101);
46  hEEpZ_ix_iy_zMap = ibooker.book2D("hEEpZ_ix_iy_zMap","", 100,1,101, 100,1,101);
47  hEEpZ_ix_iy_dxMap = ibooker.book2D("hEEpZ_ix_iy_dxMap","", 100,1,101, 100,1,101);
48  hEEpZ_ix_iy_dyMap = ibooker.book2D("hEEpZ_ix_iy_dyMap","", 100,1,101, 100,1,101);
49  // ECAL -endcap
50  hEEmZ_ix_iy_irMap = ibooker.book2D("hEEmZ_ix_iy_irMap","", 100,1,101, 100,1,101);
51  hEEmZ_ix_iy_xMap = ibooker.book2D("hEEmZ_ix_iy_xMap","", 100,1,101, 100,1,101);
52  hEEmZ_ix_iy_yMap = ibooker.book2D("hEEmZ_ix_iy_yMap","", 100,1,101, 100,1,101);
53  hEEmZ_ix_iy_zMap = ibooker.book2D("hEEmZ_ix_iy_zMap","", 100,1,101, 100,1,101);
54  hEEmZ_ix_iy_dxMap = ibooker.book2D("hEEmZ_ix_iy_dxMap","", 100,1,101, 100,1,101);
55  hEEmZ_ix_iy_dyMap = ibooker.book2D("hEEmZ_ix_iy_dyMap","", 100,1,101, 100,1,101);
56 
57  // Initialize bins for geometry to -999 because z = 0 is a valid entry
58  for (int i=1; i<=100; i++)
59  for (int j=1; j<=100; j++)
60  {
67 
74  }
75 
76  for (int i=1; i<=171; i++)
77  {
80  for (int j=1; j<=360; j++)
81  {
84  }
85  }
86 
87  // Book Data Histograms
89 
90  hECAL_Nevents = ibooker.book1D("hECAL_Nevents","",1,0,1);
91 
92 
93  // Energy Histograms by logical index
94  hEEpZ_energy_ix_iy = ibooker.book2D("hEEpZ_energy_ix_iy","", 100,1,101, 100,1,101);
95  hEEmZ_energy_ix_iy = ibooker.book2D("hEEmZ_energy_ix_iy","", 100,1,101, 100,1,101);
96  hEB_energy_ieta_iphi = ibooker.book2D("hEB_energy_ieta_iphi","", 171, -85, 86, 360, 1, 361);
97 
98  hEEpZ_Minenergy_ix_iy = ibooker.book2D("hEEpZ_Minenergy_ix_iy","", 100,1,101, 100,1,101);
99  hEEmZ_Minenergy_ix_iy = ibooker.book2D("hEEmZ_Minenergy_ix_iy","", 100,1,101, 100,1,101);
100  hEB_Minenergy_ieta_iphi = ibooker.book2D("hEB_Minenergy_ieta_iphi","", 171, -85, 86, 360, 1, 361);
101 
102  hEEpZ_Maxenergy_ix_iy = ibooker.book2D("hEEpZ_Maxenergy_ix_iy","", 100,1,101, 100,1,101);
103  hEEmZ_Maxenergy_ix_iy = ibooker.book2D("hEEmZ_Maxenergy_ix_iy","", 100,1,101, 100,1,101);
104  hEB_Maxenergy_ieta_iphi = ibooker.book2D("hEB_Maxenergy_ieta_iphi","", 171, -85, 86, 360, 1, 361);
105 
106  // need to initialize those
107  for (int i=1; i<=171; i++)
108  for (int j=1; j<=360; j++)
109  {
112  }
113  for (int i=1; i<=100; i++)
114  for (int j=1; j<=100; j++)
115  {
120  }
121 
122 
123  // Occupancy Histograms by logical index
124  hEEpZ_Occ_ix_iy = ibooker.book2D("hEEpZ_Occ_ix_iy","", 100,1,101, 100,1,101);
125  hEEmZ_Occ_ix_iy = ibooker.book2D("hEEmZ_Occ_ix_iy","", 100,1,101, 100,1,101);
126  hEB_Occ_ieta_iphi = ibooker.book2D("hEB_Occ_ieta_iphi","",171, -85, 86, 360, 1, 361);
127 
128  // Integrated Histograms
129  if(finebinning_)
130  {
131  hEEpZ_energyvsir = ibooker.book2D("hEEpZ_energyvsir","", 100,1,101, 20110,-10,201);
132  hEEmZ_energyvsir = ibooker.book2D("hEEmZ_energyvsir","", 100,1,101, 20110,-10,201);
133  hEB_energyvsieta = ibooker.book2D("hEB_energyvsieta","", 171, -85, 86, 20110, -10, 201);
134 
135  hEEpZ_Maxenergyvsir = ibooker.book2D("hEEpZ_Maxenergyvsir","", 100,1,101, 20110,-10,201);
136  hEEmZ_Maxenergyvsir = ibooker.book2D("hEEmZ_Maxenergyvsir","", 100,1,101, 20110,-10,201);
137  hEB_Maxenergyvsieta = ibooker.book2D("hEB_Maxenergyvsieta","", 171, -85, 86, 20110, -10, 201);
138 
139  hEEpZ_Minenergyvsir = ibooker.book2D("hEEpZ_Minenergyvsir","", 100,1,101, 20110,-10,201);
140  hEEmZ_Minenergyvsir = ibooker.book2D("hEEmZ_Minenergyvsir","", 100,1,101, 20110,-10,201);
141  hEB_Minenergyvsieta = ibooker.book2D("hEB_Minenergyvsieta","", 171, -85, 86, 20110, -10, 201);
142 
143  hEEpZ_SETvsir = ibooker.book2D("hEEpZ_SETvsir","", 50,1,51, 20010,0,201);
144  hEEmZ_SETvsir = ibooker.book2D("hEEmZ_SETvsir","", 50,1,51, 20010,0,201);
145  hEB_SETvsieta = ibooker.book2D("hEB_SETvsieta","", 171, -85, 86, 20010, 0, 201);
146 
147  hEEpZ_METvsir = ibooker.book2D("hEEpZ_METvsir","", 50,1,51, 20010,0,201);
148  hEEmZ_METvsir = ibooker.book2D("hEEmZ_METvsir","", 50,1,51, 20010,0,201);
149  hEB_METvsieta = ibooker.book2D("hEB_METvsieta","", 171, -85, 86, 20010, 0, 201);
150 
151  hEEpZ_METPhivsir = ibooker.book2D("hEEpZ_METPhivsir","", 50,1,51, 80,-4,4);
152  hEEmZ_METPhivsir = ibooker.book2D("hEEmZ_METPhivsir","", 50,1,51, 80,-4,4);
153  hEB_METPhivsieta = ibooker.book2D("hEB_METPhivsieta","", 171, -85, 86, 80,-4,4);
154 
155  hEEpZ_MExvsir = ibooker.book2D("hEEpZ_MExvsir","", 50,1,51, 10010,-50,51);
156  hEEmZ_MExvsir = ibooker.book2D("hEEmZ_MExvsir","", 50,1,51, 10010,-50,51);
157  hEB_MExvsieta = ibooker.book2D("hEB_MExvsieta","", 171, -85, 86, 10010,-50,51);
158 
159  hEEpZ_MEyvsir = ibooker.book2D("hEEpZ_MEyvsir","", 50,1,51, 10010,-50,51);
160  hEEmZ_MEyvsir = ibooker.book2D("hEEmZ_MEyvsir","", 50,1,51, 10010,-50,51);
161  hEB_MEyvsieta = ibooker.book2D("hEB_MEyvsieta","", 171, -85, 86, 10010,-50,51);
162 
163  hEEpZ_Occvsir = ibooker.book2D("hEEpZ_Occvsir","", 50,1,51, 1000,0,1000);
164  hEEmZ_Occvsir = ibooker.book2D("hEEmZ_Occvsir","", 50,1,51, 1000,0,1000);
165  hEB_Occvsieta = ibooker.book2D("hEB_Occvsieta","", 171, -85, 86, 400,0,400);
166  }
167  else
168  {
169  hEEpZ_energyvsir = ibooker.book2D("hEEpZ_energyvsir","", 100,1,101, 510,-10,100);
170  hEEmZ_energyvsir = ibooker.book2D("hEEmZ_energyvsir","", 100,1,101, 510,-10,100);
171  hEB_energyvsieta = ibooker.book2D("hEB_energyvsieta","", 171, -85, 86, 510, -10, 100);
172 
173  hEEpZ_Maxenergyvsir = ibooker.book2D("hEEpZ_Maxenergyvsir","", 100,1,101, 510,-10,100);
174  hEEmZ_Maxenergyvsir = ibooker.book2D("hEEmZ_Maxenergyvsir","", 100,1,101, 510,-10,100);
175  hEB_Maxenergyvsieta = ibooker.book2D("hEB_Maxenergyvsieta","", 171, -85, 86, 510, -10, 100);
176 
177  hEEpZ_Minenergyvsir = ibooker.book2D("hEEpZ_Minenergyvsir","", 100,1,101, 510,-10,100);
178  hEEmZ_Minenergyvsir = ibooker.book2D("hEEmZ_Minenergyvsir","", 100,1,101, 510,-10,100);
179  hEB_Minenergyvsieta = ibooker.book2D("hEB_Minenergyvsieta","", 171, -85, 86, 510, -10, 100);
180 
181  hEEpZ_SETvsir = ibooker.book2D("hEEpZ_SETvsir","", 50,1,51, 510,0,100);
182  hEEmZ_SETvsir = ibooker.book2D("hEEmZ_SETvsir","", 50,1,51, 510,0,100);
183  hEB_SETvsieta = ibooker.book2D("hEB_SETvsieta","", 171, -85, 86, 510, 0, 100);
184 
185  hEEpZ_METvsir = ibooker.book2D("hEEpZ_METvsir","", 50,1,51, 510,0,100);
186  hEEmZ_METvsir = ibooker.book2D("hEEmZ_METvsir","", 50,1,51, 510,0,100);
187  hEB_METvsieta = ibooker.book2D("hEB_METvsieta","", 171, -85, 86, 510, 0, 100);
188 
189  hEEpZ_METPhivsir = ibooker.book2D("hEEpZ_METPhivsir","", 50,1,51, 80,-4,4);
190  hEEmZ_METPhivsir = ibooker.book2D("hEEmZ_METPhivsir","", 50,1,51, 80,-4,4);
191  hEB_METPhivsieta = ibooker.book2D("hEB_METPhivsieta","", 171, -85, 86, 80,-4,4);
192 
193  hEEpZ_MExvsir = ibooker.book2D("hEEpZ_MExvsir","", 50,1,51, 510,-50,51);
194  hEEmZ_MExvsir = ibooker.book2D("hEEmZ_MExvsir","", 50,1,51, 510,-50,51);
195  hEB_MExvsieta = ibooker.book2D("hEB_MExvsieta","", 171, -85, 86, 510,-50,51);
196 
197  hEEpZ_MEyvsir = ibooker.book2D("hEEpZ_MEyvsir","", 50,1,51, 510,-50,51);
198  hEEmZ_MEyvsir = ibooker.book2D("hEEmZ_MEyvsir","", 50,1,51, 510,-50,51);
199  hEB_MEyvsieta = ibooker.book2D("hEB_MEyvsieta","", 171, -85, 86, 510,-50,51);
200 
201  hEEpZ_Occvsir = ibooker.book2D("hEEpZ_Occvsir","", 50,1,51, 1000,0,1000);
202  hEEmZ_Occvsir = ibooker.book2D("hEEmZ_Occvsir","", 50,1,51, 1000,0,1000);
203  hEB_Occvsieta = ibooker.book2D("hEB_Occvsieta","", 171, -85, 86, 400,0,400);
204 
205  }
206 }
207 
209 {
210  // Fill geometry histograms
211  using namespace edm;
212  //int b=0;
214  iSetup.get<CaloGeometryRecord>().get(pG);
215  const CaloGeometry cG = *pG;
216 
217  //----Fill Ecal Barrel----//
219  int n=0;
220  std::vector<DetId> EBids=EBgeom->getValidDetIds(DetId::Ecal, 1);
221  for (std::vector<DetId>::iterator i=EBids.begin(); i!=EBids.end(); i++) {
222  n++;
223  auto cell=EBgeom->getGeometry(*i);
224  //GlobalPoint p = cell->getPosition();
225 
226  EBDetId EcalID(i->rawId());
227 
228  int Crystal_ieta = EcalID.ieta();
229  int Crystal_iphi = EcalID.iphi();
230  double Crystal_eta = cell->getPosition().eta();
231  double Crystal_phi = cell->getPosition().phi();
232  hEB_ieta_iphi_etaMap->setBinContent(Crystal_ieta+86, Crystal_iphi, Crystal_eta);
233  hEB_ieta_iphi_phiMap->setBinContent(Crystal_ieta+86, Crystal_iphi, (Crystal_phi*180/M_PI) );
234 
235  DEBUG( " Crystal " << n );
236  DEBUG( " ieta, iphi = " << Crystal_ieta << ", " << Crystal_iphi);
237  DEBUG( " eta, phi = " << cell->getPosition().eta() << ", " << cell->getPosition().phi());
238  DEBUG( " " );
239 
240  }
241  //----Fill Ecal Endcap----------//
243  n=0;
244  std::vector<DetId> EEids=EEgeom->getValidDetIds(DetId::Ecal, 2);
245  for (std::vector<DetId>::iterator i=EEids.begin(); i!=EEids.end(); i++) {
246  n++;
247  auto cell=EEgeom->getGeometry(*i);
248  //GlobalPoint p = cell->getPosition();
249  EEDetId EcalID(i->rawId());
250  int Crystal_zside = EcalID.zside();
251  int Crystal_ix = EcalID.ix();
252  int Crystal_iy = EcalID.iy();
253  Float_t ix_ = Crystal_ix-50.5;
254  Float_t iy_ = Crystal_iy-50.5;
255  Int_t ir = (Int_t)sqrt(ix_*ix_ + iy_*iy_);
256 
257  //double Crystal_eta = cell->getPosition().eta();
258  //double Crystal_phi = cell->getPosition().phi();
259  double Crystal_x = cell->getPosition().x();
260  double Crystal_y = cell->getPosition().y();
261  double Crystal_z = cell->getPosition().z();
262  // ECAL -endcap
263  if (Crystal_zside == -1)
264  {
265  hEEmZ_ix_iy_irMap->setBinContent(Crystal_ix, Crystal_iy, ir);
266  hEEmZ_ix_iy_xMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_x);
267  hEEmZ_ix_iy_yMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_y);
268  hEEmZ_ix_iy_zMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_z);
269  }
270  // ECAL +endcap
271  if (Crystal_zside == 1)
272  {
273  hEEpZ_ix_iy_irMap->setBinContent(Crystal_ix, Crystal_iy, ir);
274  hEEpZ_ix_iy_xMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_x);
275  hEEpZ_ix_iy_yMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_y);
276  hEEpZ_ix_iy_zMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_z);
277  }
278 
279  DEBUG( " Crystal " << n );
280  DEBUG( " side = " << Crystal_zside );
281  DEBUG(" ix, iy = " << Crystal_ix << ", " << Crystal_iy);
282  DEBUG(" x, y = " << Crystal_x << ", " << Crystal_y);;
283  DEBUG( " " );
284 
285  }
286 
287  //-------Set the cell size for each (ieta, iphi) bin-------//
288  double currentLowEdge_eta = 0;
289  //double currentHighEdge_eta = 0;
290  for (int ieta=1; ieta<=85 ; ieta++)
291  {
292  int ieta_ = 86 + ieta;
293 
294  double eta = hEB_ieta_iphi_etaMap->getBinContent(ieta_, 1);
295  double etam1 = -999;
296 
297  if (ieta==1)
298  etam1 = hEB_ieta_iphi_etaMap->getBinContent(85, 1);
299  else
300  etam1 = hEB_ieta_iphi_etaMap->getBinContent(ieta_ - 1, 1);
301 
302  //double phi = hEB_ieta_iphi_phiMap->getBinContent(ieta_, 1);
303  double deta = fabs( eta - etam1 );
304  double dphi = fabs( hEB_ieta_iphi_phiMap->getBinContent(ieta_, 1) - hEB_ieta_iphi_phiMap->getBinContent(ieta_, 2) );
305 
306  currentLowEdge_eta += deta;
307  hEB_ieta_detaMap->setBinContent(ieta_, deta); // positive rings
308  hEB_ieta_dphiMap->setBinContent(ieta_, dphi); // positive rings
309  hEB_ieta_detaMap->setBinContent(86-ieta, deta); // negative rings
310  hEB_ieta_dphiMap->setBinContent(86-ieta, dphi); // negative rings
311  }
312 }
313 
314 
315 
317 {
318  CurrentEvent++;
319  DEBUG( "Event: " << CurrentEvent);
320  WriteECALRecHits( iEvent, iSetup );
321  hECAL_Nevents->Fill(0.5);
322 }
323 
325 {
328  iEvent.getByToken( EBRecHitsLabel_, EBRecHits );
329  iEvent.getByToken( EERecHitsLabel_, EERecHits );
330  DEBUG( "Got ECALRecHits");
331 
333  iSetup.get<CaloGeometryRecord>().get(pG);
334  const CaloGeometry cG = *pG;
335  //const CaloSubdetectorGeometry* EBgeom=cG.getSubdetectorGeometry(DetId::Ecal,1);
336  //const CaloSubdetectorGeometry* EEgeom=cG.getSubdetectorGeometry(DetId::Ecal,2);
337  DEBUG( "Got Geometry");
338 
339  TLorentzVector vEBMET_EtaRing[171];
340  int EBActiveRing[171];
341  int EBNActiveCells[171];
342  double EBSET_EtaRing[171];
343  double EBMaxEnergy_EtaRing[171];
344  double EBMinEnergy_EtaRing[171];
345  double EBenergy_EtaRing[171];
346 
347  for (int i=0; i<171; i++)
348  {
349  EBActiveRing[i] = 0;
350  EBNActiveCells[i] = 0;
351  EBSET_EtaRing[i] = 0.0;
352  EBMaxEnergy_EtaRing[i] = -999;
353  EBMinEnergy_EtaRing[i] = 14E3;
354  EBenergy_EtaRing[i] = 0.0;
355  }
356 
357  edm::LogInfo("OutputInfo") << "Looping over EB" << std::endl;
358 
360  //int nEBrechit = 0;
361 
362  for (ebrechit = EBRecHits->begin(); ebrechit != EBRecHits->end(); ebrechit++) {
363 
364  EBDetId det = ebrechit->id();
365  double Energy = ebrechit->energy();
366  Int_t ieta = det.ieta();
367  Int_t iphi = det.iphi();
368  int EtaRing = 85 + ieta; // this counts from 0
369  double eta = hEB_ieta_iphi_etaMap->getBinContent(EtaRing+1,iphi);
370  double phi = hEB_ieta_iphi_phiMap->getBinContent(EtaRing+1,iphi);
371  double theta = 2*TMath::ATan(exp(-1*eta));
372  double ET = Energy*TMath::Sin(theta);
373  TLorentzVector v_;
374 
375  if (Energy>EBMaxEnergy_EtaRing[EtaRing]) EBMaxEnergy_EtaRing[EtaRing] = Energy;
376  if (Energy<EBMinEnergy_EtaRing[EtaRing]) EBMinEnergy_EtaRing[EtaRing] = Energy;
377 
378  if (Energy>0)
379  {
380  EBActiveRing[EtaRing] = 1;
381  EBNActiveCells[EtaRing]++;
382  EBSET_EtaRing[EtaRing]+=ET;
383  v_.SetPtEtaPhiE(ET, 0, phi, ET);
384  vEBMET_EtaRing[EtaRing]-=v_;
385  EBenergy_EtaRing[EtaRing]+=Energy;
386  hEB_Occ_ieta_iphi->Fill(ieta, iphi);
387 
388  }
389 
390  hEB_energy_ieta_iphi->Fill(ieta, iphi, Energy);
391  if (Energy>hEB_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi))
392  hEB_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
393  if (Energy<hEB_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi))
394  hEB_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
395 
396 
397  } // loop over EB
398 
399  for (int iEtaRing = 0; iEtaRing < 171; iEtaRing++)
400  {
401  hEB_Minenergyvsieta->Fill(iEtaRing-85, EBMinEnergy_EtaRing[iEtaRing]);
402  hEB_Maxenergyvsieta->Fill(iEtaRing-85, EBMaxEnergy_EtaRing[iEtaRing]);
403 
404  if (EBActiveRing[iEtaRing])
405  {
406  hEB_METvsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Pt());
407  hEB_METPhivsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Phi());
408  hEB_MExvsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Px());
409  hEB_MEyvsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Py());
410  hEB_SETvsieta->Fill(iEtaRing-85, EBSET_EtaRing[iEtaRing]);
411  hEB_Occvsieta->Fill(iEtaRing-85, EBNActiveCells[iEtaRing]);
412  hEB_energyvsieta->Fill(iEtaRing-85, EBenergy_EtaRing[iEtaRing]);
413  }
414  }
415 
416 
417  TLorentzVector vEEpZMET_EtaRing[101];
418  int EEpZActiveRing[101];
419  int EEpZNActiveCells[101];
420  double EEpZSET_EtaRing[101];
421  double EEpZMaxEnergy_EtaRing[101];
422  double EEpZMinEnergy_EtaRing[101];
423 
424  TLorentzVector vEEmZMET_EtaRing[101];
425  int EEmZActiveRing[101];
426  int EEmZNActiveCells[101];
427  double EEmZSET_EtaRing[101];
428  double EEmZMaxEnergy_EtaRing[101];
429  double EEmZMinEnergy_EtaRing[101];
430 
431  for (int i=0;i<101; i++)
432  {
433  EEpZActiveRing[i] = 0;
434  EEpZNActiveCells[i] = 0;
435  EEpZSET_EtaRing[i] = 0.0;
436  EEpZMaxEnergy_EtaRing[i] = -999;
437  EEpZMinEnergy_EtaRing[i] = 14E3;
438 
439  EEmZActiveRing[i] = 0;
440  EEmZNActiveCells[i] = 0;
441  EEmZSET_EtaRing[i] = 0.0;
442  EEmZMaxEnergy_EtaRing[i] = -999;
443  EEmZMinEnergy_EtaRing[i] = 14E3;
444  }
445 
446  edm::LogInfo("OutputInfo") << "Looping over EE" << std::endl;
448  //int nEErechit = 0;
449  for (eerechit = EERecHits->begin(); eerechit != EERecHits->end(); eerechit++) {
450 
451  EEDetId det = eerechit->id();
452  double Energy = eerechit->energy();
453  Int_t ix = det.ix();
454  Int_t iy = det.iy();
455  //Float_t ix_ = (Float_t)-999;
456  //Float_t iy_ = (Float_t)-999;
457  Int_t ir = -999;
458  // edm::LogInfo("OutputInfo") << ix << " " << iy << " " << ix_ << " " << iy_ << " " << ir << std::endl;
459 
460  double x = -999;
461  double y = -999;
462  double z = -999;
463  double theta = -999;
464  double phi = -999;
465 
466  int Crystal_zside = det.zside();
467 
468  if (Crystal_zside == -1)
469  {
470  ir = (Int_t)hEEmZ_ix_iy_irMap->getBinContent(ix,iy);
471  x = hEEmZ_ix_iy_xMap->getBinContent(ix,iy);
472  y = hEEmZ_ix_iy_yMap->getBinContent(ix,iy);
473  z = hEEmZ_ix_iy_zMap->getBinContent(ix,iy);
474  }
475  if (Crystal_zside == 1)
476  {
477  ir = (Int_t)hEEpZ_ix_iy_irMap->getBinContent(ix,iy);
478  x = hEEpZ_ix_iy_xMap->getBinContent(ix,iy);
479  y = hEEpZ_ix_iy_yMap->getBinContent(ix,iy);
480  z = hEEpZ_ix_iy_zMap->getBinContent(ix,iy);
481  }
482  TVector3 pos_vector(x,y,z);
483  phi = pos_vector.Phi();
484  theta = pos_vector.Theta();
485  double ET = Energy*TMath::Sin(theta);
486  TLorentzVector v_;
487 
488 
489  if (Crystal_zside == -1)
490  {
491  if (Energy>0)
492  {
493  EEmZActiveRing[ir] = 1;
494  EEmZNActiveCells[ir]++;
495  EEmZSET_EtaRing[ir]+=ET;
496  v_.SetPtEtaPhiE(ET,0,phi,ET);
497  vEEmZMET_EtaRing[ir]-=v_;
498  hEEmZ_Occ_ix_iy->Fill(ix, iy);
499  }
500  hEEmZ_energyvsir->Fill(ir, Energy);
501  hEEmZ_energy_ix_iy->Fill(ix, iy, Energy);
502 
503  if (Energy>EEmZMaxEnergy_EtaRing[ir]) EEmZMaxEnergy_EtaRing[ir] = Energy;
504  if (Energy<EEmZMinEnergy_EtaRing[ir]) EEmZMinEnergy_EtaRing[ir] = Energy;
505 
506  if (Energy>hEEmZ_Maxenergy_ix_iy->getBinContent(ix,iy))
507  hEEmZ_Maxenergy_ix_iy->setBinContent(ix,iy, Energy);
508  if (Energy<hEEmZ_Minenergy_ix_iy->getBinContent(ix,iy))
509  hEEmZ_Minenergy_ix_iy->setBinContent(ix,iy, Energy);
510  }
511  if (Crystal_zside == 1)
512  {
513  if (Energy>0)
514  {
515  EEpZActiveRing[ir] = 1;
516  EEpZNActiveCells[ir]++;
517  EEpZSET_EtaRing[ir]+=ET;
518  v_.SetPtEtaPhiE(ET,0,phi,ET);
519  vEEpZMET_EtaRing[ir]-=v_;
520  hEEpZ_Occ_ix_iy->Fill(ix, iy);
521  }
522  hEEpZ_energyvsir->Fill(ir, Energy);
523  hEEpZ_energy_ix_iy->Fill(ix, iy, Energy);
524 
525  if (Energy>EEpZMaxEnergy_EtaRing[ir]) EEpZMaxEnergy_EtaRing[ir] = Energy;
526  if (Energy<EEpZMinEnergy_EtaRing[ir]) EEpZMinEnergy_EtaRing[ir] = Energy;
527  if (Energy>hEEpZ_Maxenergy_ix_iy->getBinContent(ix,iy))
528  hEEpZ_Maxenergy_ix_iy->setBinContent(ix,iy, Energy);
529  if (Energy<hEEpZ_Minenergy_ix_iy->getBinContent(ix,iy))
530  hEEpZ_Minenergy_ix_iy->setBinContent(ix,iy, Energy);
531  }
532  } // loop over EE
533  edm::LogInfo("OutputInfo") << "Done Looping over EE" << std::endl;
534  for (int iEtaRing = 0; iEtaRing<101; iEtaRing++)
535  {
536  hEEpZ_Maxenergyvsir->Fill(iEtaRing, EEpZMaxEnergy_EtaRing[iEtaRing]);
537  hEEpZ_Minenergyvsir->Fill(iEtaRing, EEpZMinEnergy_EtaRing[iEtaRing]);
538  hEEmZ_Maxenergyvsir->Fill(iEtaRing, EEmZMaxEnergy_EtaRing[iEtaRing]);
539  hEEmZ_Minenergyvsir->Fill(iEtaRing, EEmZMinEnergy_EtaRing[iEtaRing]);
540 
541  if (EEpZActiveRing[iEtaRing])
542  {
543  hEEpZ_METvsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Pt());
544  hEEpZ_METPhivsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Phi());
545  hEEpZ_MExvsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Px());
546  hEEpZ_MEyvsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Py());
547  hEEpZ_SETvsir->Fill(iEtaRing, EEpZSET_EtaRing[iEtaRing]);
548  hEEpZ_Occvsir->Fill(iEtaRing, EEpZNActiveCells[iEtaRing]);
549  }
550 
551  if (EEmZActiveRing[iEtaRing])
552  {
553  hEEmZ_METvsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Pt());
554  hEEmZ_METPhivsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Phi());
555  hEEmZ_MExvsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Px());
556  hEEmZ_MEyvsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Py());
557  hEEmZ_SETvsir->Fill(iEtaRing, EEmZSET_EtaRing[iEtaRing]);
558  hEEmZ_Occvsir->Fill(iEtaRing, EEmZNActiveCells[iEtaRing]);
559  }
560  }
561  edm::LogInfo("OutputInfo") << "Done ..." << std::endl;
562 } // loop over RecHits
563 
564 
565 
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:44
T getUntrackedParameter(std::string const &, T const &) const
#define DEBUG(X)
MonitorElement * hEB_Maxenergyvsieta
MonitorElement * hEB_MExvsieta
MonitorElement * hEEmZ_ix_iy_dxMap
void setBinContent(int binx, double content)
set content of bin (1-D)
int ix() const
Definition: EEDetId.h:76
MonitorElement * hEEmZ_Occvsir
MonitorElement * hEEpZ_METPhivsir
MonitorElement * hEEmZ_METvsir
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:519
MonitorElement * hEEpZ_METvsir
MonitorElement * hEB_Maxenergy_ieta_iphi
MonitorElement * hEB_METvsieta
MonitorElement * hEEpZ_ix_iy_irMap
MonitorElement * hEB_Occvsieta
std::vector< EcalRecHit >::const_iterator const_iterator
Geom::Theta< T > theta() const
MonitorElement * hEB_energyvsieta
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
MonitorElement * hEEpZ_ix_iy_zMap
MonitorElement * hEEpZ_Maxenergy_ix_iy
MonitorElement * hEB_METPhivsieta
virtual const std::vector< DetId > & getValidDetIds(DetId::Detector det=DetId::Detector(0), int subdet=0) const
Get a list of valid detector ids (for the given subdetector)
MonitorElement * hEEmZ_Minenergyvsir
void Fill(long long x)
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
int iEvent
Definition: GenABIO.cc:230
MonitorElement * hEEmZ_energy_ix_iy
MonitorElement * hEEmZ_energyvsir
MonitorElement * hEEmZ_ix_iy_dyMap
MonitorElement * hEEmZ_Occ_ix_iy
MonitorElement * hECAL_Nevents
MonitorElement * hEB_ieta_dphiMap
T sqrt(T t)
Definition: SSEVec.h:18
MonitorElement * hEEmZ_Maxenergyvsir
MonitorElement * hEEmZ_SETvsir
MonitorElement * hEB_Occ_ieta_iphi
MonitorElement * hEEpZ_Minenergyvsir
int zside() const
Definition: EEDetId.h:70
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:118
MonitorElement * hEEpZ_Minenergy_ix_iy
MonitorElement * hEB_Minenergy_ieta_iphi
MonitorElement * hEEmZ_METPhivsir
MonitorElement * hEEmZ_MEyvsir
MonitorElement * hEEpZ_energyvsir
void FillGeometry(const edm::EventSetup &)
MonitorElement * hEEpZ_ix_iy_dyMap
int iy() const
Definition: EEDetId.h:82
MonitorElement * hEB_ieta_detaMap
MonitorElement * hEEmZ_Minenergy_ix_iy
int ieta() const
get the crystal ieta
Definition: EBDetId.h:51
MonitorElement * hEB_energy_ieta_iphi
MonitorElement * hEB_MEyvsieta
#define M_PI
MonitorElement * hEB_Minenergyvsieta
const_iterator end() const
ECALRecHitAnalyzer(const edm::ParameterSet &)
MonitorElement * hEEpZ_ix_iy_yMap
MonitorElement * hEEmZ_ix_iy_zMap
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:279
MonitorElement * hEEmZ_Maxenergy_ix_iy
virtual void dqmbeginRun(const edm::Run &, const edm::EventSetup &)
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:136
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
const T & get() const
Definition: EventSetup.h:59
MonitorElement * hEB_ieta_iphi_etaMap
MonitorElement * hEEpZ_SETvsir
MonitorElement * hEEpZ_Maxenergyvsir
double getBinContent(int binx) const
get content of bin (1-D)
HLT enums.
MonitorElement * hEEmZ_MExvsir
MonitorElement * hEEmZ_ix_iy_yMap
MonitorElement * hEEmZ_ix_iy_xMap
MonitorElement * hEEpZ_MEyvsir
void WriteECALRecHits(const edm::Event &, const edm::EventSetup &)
MonitorElement * hEB_SETvsieta
MonitorElement * hEEpZ_MExvsir
MonitorElement * hEEpZ_Occvsir
MonitorElement * hEEpZ_energy_ix_iy
MonitorElement * hEEpZ_Occ_ix_iy
#define ET
MonitorElement * hEEmZ_ix_iy_irMap
edm::EDGetTokenT< EERecHitCollection > EERecHitsLabel_
MonitorElement * hEB_ieta_iphi_phiMap
MonitorElement * hEEpZ_ix_iy_xMap
edm::EDGetTokenT< EBRecHitCollection > EBRecHitsLabel_
const_iterator begin() const
Definition: Run.h:43
void analyze(const edm::Event &, const edm::EventSetup &) override
MonitorElement * hEEpZ_ix_iy_dxMap