CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Analyzer_minbias.cc
Go to the documentation of this file.
1 // system include files
2 #include <memory>
3 #include <string>
4 #include <iostream>
5 
6 // user include files
7 #include "Geometry/Records/interface/IdealGeometryRecord.h"
8 #include "Calibration/HcalCalibAlgos/interface/Analyzer_minbias.h"
9 #include "DataFormats/Provenance/interface/Provenance.h"
10 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
11 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
12 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
13 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
14 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
15 #include "DataFormats/L1GlobalTrigger/interface/L1GtfeWord.h"
16 #include "FWCore/Framework/interface/Run.h"
17 
18 #include "TFile.h"
19 #include "TH1.h"
20 #include "TH2.h"
21 #include <fstream>
22 #include <sstream>
23 
24 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
25 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
26 
27 using namespace std;
28 using namespace reco;
29 //
30 // constructors and destructor
31 //
32 namespace cms{
33 Analyzer_minbias::Analyzer_minbias(const edm::ParameterSet& iConfig)
34 {
35  // get name of output file with histogramms
36  fOutputFileName = iConfig.getUntrackedParameter<string>("HistOutFile");
37  // get names of modules, producing object collections
38 
39  hbherecoMB = iConfig.getParameter<edm::InputTag>("hbheInputMB");
40  horecoMB = iConfig.getParameter<edm::InputTag>("hoInputMB");
41  hfrecoMB = iConfig.getParameter<edm::InputTag>("hfInputMB");
42 
43  hbherecoNoise = iConfig.getParameter<edm::InputTag>("hbheInputNoise");
44  horecoNoise = iConfig.getParameter<edm::InputTag>("hoInputNoise");
45  hfrecoNoise = iConfig.getParameter<edm::InputTag>("hfInputNoise");
46 
47  theRecalib = iConfig.getParameter<bool>("Recalib");
48 
49 //
50 //
51  for(int i=0; i<73; i++)
52  {
53  for(int j=0; j<43; j++)
54  {
55  noise_min[i][j] = 0.;
56  noise_pl[i][j] = 0.;
57  }
58  }
59 //
60 //
61 
62 }
63 
64 Analyzer_minbias::~Analyzer_minbias()
65 {
66 
67  // do anything here that needs to be done at desctruction time
68  // (e.g. close files, deallocate resources etc.)
69 
70 }
71 
72 void Analyzer_minbias::beginRun( const edm::Run& r, const edm::EventSetup& iSetup)
73 {
74  nevent_run = 0;
75 }
76 void Analyzer_minbias::endRun( const edm::Run& r, const edm::EventSetup& iSetup)
77 {
78  std::cout<<" Runnumber "<<r.run()<<" Nevents "<<nevent_run<<std::endl;
79 }
80 
81 void Analyzer_minbias::beginJob()
82 {
83 
84  hOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
85 
86  myTree = new TTree("RecJet","RecJet Tree");
87  myTree->Branch("mydet", &mydet, "mydet/I");
88  myTree->Branch("mysubd", &mysubd, "mysubd/I");
89  myTree->Branch("depth", &depth, "depth/I");
90  myTree->Branch("ieta", &ieta, "ieta/I");
91  myTree->Branch("iphi", &iphi, "iphi/I");
92  myTree->Branch("eta", &eta, "eta/F");
93  myTree->Branch("phi", &phi, "phi/F");
94 
95  myTree->Branch("mom0_MB", &mom0_MB, "mom0_MB/F");
96  myTree->Branch("mom1_MB", &mom1_MB, "mom1_MB/F");
97  myTree->Branch("mom2_MB", &mom2_MB, "mom2_MB/F");
98  myTree->Branch("mom4_MB", &mom4_MB, "mom4_MB/F");
99 
100  myTree->Branch("mom0_Noise", &mom0_Noise, "mom0_Noise/F");
101  myTree->Branch("mom1_Noise", &mom1_Noise, "mom1_Noise/F");
102  myTree->Branch("mom2_Noise", &mom2_Noise, "mom2_Noise/F");
103  myTree->Branch("mom4_Noise", &mom4_Noise, "mom4_Noise/F");
104 
105  myTree->Branch("mom0_Diff", &mom0_Diff, "mom0_Diff/F");
106  myTree->Branch("mom1_Diff", &mom1_Diff, "mom1_Diff/F");
107  myTree->Branch("mom2_Diff", &mom2_Diff, "mom2_Diff/F");
108 
109  myTree->Branch("occup", &occup, "occup/F");
110 
111  std::cout<<" Before ordering Histos "<<std::endl;
112 
113  char str0[15];
114  char str1[15];
115 
116  char str10[15];
117  char str11[15];
118 
119  int k=0;
120  nevent = 0;
121 // Size of collections
122 
123  hHBHEsize_vs_run = new TH2F("hHBHEsize_vs_run","hHBHEsize_vs_run",500,111500.,112000.,6101,-100.5,6000.5);
124  hHFsize_vs_run = new TH2F("hHFsize_vs_run","hHFsize_vs_run",500,111500.,112000.,6101,-100.5,6000.5);
125 
126  for(int i=1;i<73;i++){
127  for(int j=1;j<43;j++){
128 
129  meannoise_pl[i][j] = 0.;
130  meannoise_min[i][j] = 0.;
131 
132 // for(int l=1;l<5;l++){
133  k = i*1000+j;
134  sprintf(str0,"mpl%d",k);
135  sprintf(str1,"mmin%d",k);
136 
137  sprintf(str10,"vpl%d",k);
138  sprintf(str11,"vmin%d",k);
139 // cout<<" "<<i<<" "<<j<<endl;
140  if( j < 30 )
141  {
142 // first order moment
143  hCalo1[i][j] = new TH1F(str0, "h0", 320, -10., 10.);
144  hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
145 
146 // second order moment
147  hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 20.);
148  hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 20.);
149  }
150  else
151  {
152 // HF
153 // first order moment
154 // cout<<" "<<i<<" "<<j<<" "<<k<<endl;
155  if(j < 40)
156  {
157  hCalo1[i][j] = new TH1F(str0, "h0", 320, -10., 10.);
158  hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
159 //
160 // second order moment
161  hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 40.);
162  hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 40.);
163  }
164  else
165  {
166  hCalo1[i][j] = new TH1F(str0,"h0" , 320, -10., 10.);
167  hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
168 
169 // second order moment
170  hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 120.);
171  hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 120.);
172 
173  }
174  } // HE/HF boundary
175 // } // l
176  } // j
177  } // i
178 
179 
180  std::cout<<" After ordering Histos "<<std::endl;
181 
182  std::string ccc = "noise_0.dat";
183 
184  myout_hcal = new ofstream(ccc.c_str());
185  if(!myout_hcal) cout << " Output file not open!!! "<<endl;
186 
187 //
188  for (int i=0; i<5;i++)
189  {
190  for (int j=0; j<5;j++)
191  {
192  for (int k=0; k<73;k++)
193  {
194  for (int l=0; l<43;l++)
195  {
196  theMBFillDetMapPl0[i][j][k][l] = 0.;
197  theMBFillDetMapPl1[i][j][k][l] = 0.;
198  theMBFillDetMapPl2[i][j][k][l] = 0.;
199  theMBFillDetMapPl4[i][j][k][l] = 0.;
200 
201  theMBFillDetMapMin0[i][j][k][l] = 0.;
202  theMBFillDetMapMin1[i][j][k][l] = 0.;
203  theMBFillDetMapMin2[i][j][k][l] = 0.;
204  theMBFillDetMapMin4[i][j][k][l] = 0.;
205 
206 
207  theNSFillDetMapPl0[i][j][k][l] = 0.;
208  theNSFillDetMapPl1[i][j][k][l] = 0.;
209  theNSFillDetMapPl2[i][j][k][l] = 0.;
210  theNSFillDetMapPl4[i][j][k][l] = 0.;
211 
212  theNSFillDetMapMin0[i][j][k][l] = 0.;
213  theNSFillDetMapMin1[i][j][k][l] = 0.;
214  theNSFillDetMapMin2[i][j][k][l] = 0.;
215  theNSFillDetMapMin4[i][j][k][l] = 0.;
216 
217  theDFFillDetMapPl0[i][j][k][l] = 0.;
218  theDFFillDetMapPl1[i][j][k][l] = 0.;
219  theDFFillDetMapPl2[i][j][k][l] = 0.;
220  theDFFillDetMapMin0[i][j][k][l] = 0.;
221  theDFFillDetMapMin1[i][j][k][l] = 0.;
222  theDFFillDetMapMin2[i][j][k][l] = 0.;
223  }
224  }
225  }
226  }
227 
228  return ;
229 }
230 //
231 // EndJob
232 //
233 void Analyzer_minbias::endJob()
234 {
235  int ii=0;
236 
237  for (int i=1; i<5;i++)
238  {
239  for (int j=1; j<5;j++)
240  {
241  for (int k=1; k<73;k++)
242  {
243  for (int l=1; l<43;l++)
244  {
245  if(theMBFillDetMapPl0[i][j][k][l] > 0)
246  {
247  mom0_MB = theMBFillDetMapPl0[i][j][k][l];
248  mom1_MB = theMBFillDetMapPl1[i][j][k][l];
249  mom2_MB = theMBFillDetMapPl2[i][j][k][l];
250  mom4_MB = theMBFillDetMapPl4[i][j][k][l];
251  mom0_Noise = theNSFillDetMapPl0[i][j][k][l];
252  mom1_Noise = theNSFillDetMapPl1[i][j][k][l];
253  mom2_Noise = theNSFillDetMapPl2[i][j][k][l];
254  mom4_Noise = theNSFillDetMapPl4[i][j][k][l];
255  mom0_Diff = theDFFillDetMapPl0[i][j][k][l];
256  mom1_Diff = theDFFillDetMapPl1[i][j][k][l];
257  mom2_Diff = theDFFillDetMapPl2[i][j][k][l];
258 
259  mysubd = i;
260  depth = j;
261  ieta = l;
262  iphi = k;
263  cout<<" Result Plus= "<<mysubd<<" "<<ieta<<" "<<iphi<<" mom0 "<<mom0_MB<<" mom1 "<<mom1_MB<<" mom2 "<<mom2_MB<<endl;
264  myTree->Fill();
265  ii++;
266  } // Pl > 0
267 
268 
269  if(theMBFillDetMapMin0[i][j][k][l] > 0)
270  {
271  mom0_MB = theMBFillDetMapMin0[i][j][k][l];
272  mom1_MB = theMBFillDetMapMin1[i][j][k][l];
273  mom2_MB = theMBFillDetMapMin2[i][j][k][l];
274  mom4_MB = theMBFillDetMapMin4[i][j][k][l];
275  mom0_Noise = theNSFillDetMapMin0[i][j][k][l];
276  mom1_Noise = theNSFillDetMapMin1[i][j][k][l];
277  mom2_Noise = theNSFillDetMapMin2[i][j][k][l];
278  mom4_Noise = theNSFillDetMapMin4[i][j][k][l];
279  mom0_Diff = theDFFillDetMapMin0[i][j][k][l];
280  mom1_Diff = theDFFillDetMapMin1[i][j][k][l];
281  mom2_Diff = theDFFillDetMapMin2[i][j][k][l];
282 
283 
284  mysubd = i;
285  depth = j;
286  ieta = -1*l;
287  iphi = k;
288  cout<<" Result Minus= "<<mysubd<<" "<<ieta<<" "<<iphi<<" mom0 "<<mom0_MB<<" mom1 "<<mom1_MB<<" mom2 "<<mom2_MB<<endl;
289  myTree->Fill();
290  ii++;
291 
292  } // Min>0
293  } // ieta
294  } // iphi
295  } // depth
296  } //subd
297 
298 
299 
300  cout<<" Number of cells "<<ii<<endl;
301 
302  hOutputFile->Write() ;
303  hOutputFile->cd();
304  hHBHEsize_vs_run->Write() ;
305  hHFsize_vs_run->Write() ;
306  for(int i=1;i<73;i++){
307  for(int j=1;j<43;j++){
308  hCalo1[i][j]->Write();
309  hCalo2[i][j]->Write();
310  hCalo1mom2[i][j]->Write();
311  hCalo2mom2[i][j]->Write();
312  }
313  }
314 
315 
316  myTree->Write();
317  hOutputFile->Close() ;
318 
319  return ;
320 }
321 
322 
323 //
324 // member functions
325 //
326 
327 // ------------ method called to produce the data ------------
328 void
329 Analyzer_minbias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
330 {
331 
332  std::cout<<" Start Analyzer_minbias::analyze "<<nevent<<std::endl;
333  nevent++;
334  nevent_run++;
335  using namespace edm;
336 
337  float rnnum = (float)iEvent.run();
338 /*
339  std::vector<Provenance const*> theProvenance;
340  iEvent.getAllProvenance(theProvenance);
341 
342  for( std::vector<Provenance const*>::const_iterator ip = theProvenance.begin();
343  ip != theProvenance.end(); ip++)
344  {
345  cout<<" Print all module/label names "<<(**ip).moduleName()<<" "<<(**ip).moduleLabel()<<
346  " "<<(**ip).productInstanceName()<<endl;
347  }
348 */
349 
350  edm::Handle<L1GlobalTriggerReadoutRecord> gtRecord;
351  iEvent.getByLabel("gtDigis", gtRecord);
352 
353  if (!gtRecord.isValid()) {
354 
355 // LogDebug("L1GlobalTriggerRecordProducer")
356 // << "\n\n Error: no L1GlobalTriggerReadoutRecord found with input tag "
357 // << m_l1GtReadoutRecord
358 // << "\n Returning empty L1GlobalTriggerRecord.\n\n"
359 // << std::endl;
360  cout<<" No L1 trigger record "<<endl;
361 // return;
362  }
363 
364 
365  const HcalRespCorrs* myRecalib=0;
366  if( theRecalib ) {
367 // Radek:
368  edm::ESHandle <HcalRespCorrs> recalibCorrs;
369  iSetup.get<HcalRespCorrsRcd>().get("recalibrate",recalibCorrs);
370  myRecalib = recalibCorrs.product();
371 // end
372  } // theRecalib
373 
374 // Noise part for HB HE
375 
376  double tmpNSFillDetMapPl1[5][5][73][43];
377  double tmpNSFillDetMapPl2[5][5][73][43];
378  double tmpNSFillDetMapMin1[5][5][73][43];
379  double tmpNSFillDetMapMin2[5][5][73][43];
380 
381  for (int i=0; i<5;i++)
382  {
383  for (int j=0; j<5;j++)
384  {
385  for (int k=0; k<73;k++)
386  {
387  for (int l=0; l<43;l++)
388  {
389  tmpNSFillDetMapPl1[i][j][k][l] = 0.;
390  tmpNSFillDetMapPl2[i][j][k][l] = 0.;
391  tmpNSFillDetMapMin1[i][j][k][l] = 0.;
392  tmpNSFillDetMapMin2[i][j][k][l] = 0.;
393  }
394  }
395  }
396  }
397  edm::Handle<HBHERecHitCollection> hbheNormal;
398  iEvent.getByLabel("hbhereco", hbheNormal);
399  if(!hbheNormal.isValid()){
400  cout<<" hbheNormal failed "<<endl;
401  } else {
402  cout<<" The size of the normal collection "<<hbheNormal->size()<<endl;
403  }
404 
405 
406  edm::Handle<HBHERecHitCollection> hbheNS;
407  iEvent.getByLabel(hbherecoNoise, hbheNS);
408 
409 
410  if(!hbheNS.isValid()){
411  LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
412  cout<<" No HBHE MS "<<endl;
413  return ;
414  }
415 
416 
417  const HBHERecHitCollection HithbheNS = *(hbheNS.product());
418  cout<<" HBHE NS size of collection "<<HithbheNS.size()<<endl;
419  hHBHEsize_vs_run->Fill(rnnum,(float)HithbheNS.size());
420 
421  if(HithbheNS.size()!= 5184) {
422  cout<<" HBHE problem "<<rnnum<<" "<<HithbheNS.size()<<endl;
423  return;
424  }
425  edm::Handle<HBHERecHitCollection> hbheMB;
426  iEvent.getByLabel(hbherecoMB, hbheMB);
427 
428  if(!hbheMB.isValid()){
429  LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
430  cout<<" No HBHE MB"<<endl;
431  return ;
432  }
433 
434  const HBHERecHitCollection HithbheMB = *(hbheMB.product());
435  cout<<" HBHE MB size of collection "<<HithbheMB.size()<<endl;
436  if(HithbheMB.size()!= 5184) {
437  cout<<" HBHE problem "<<rnnum<<" "<<HithbheMB.size()<<endl;
438  return;
439  }
440 
441  edm::Handle<HFRecHitCollection> hfNS;
442  iEvent.getByLabel(hfrecoNoise, hfNS);
443 
444  if(!hfNS.isValid()){
445  LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
446  cout<<" No HF NS "<<endl;
447  return ;
448  }
449 
450  const HFRecHitCollection HithfNS = *(hfNS.product());
451  cout<<" HFE NS size of collection "<<HithfNS.size()<<endl;
452  hHFsize_vs_run->Fill(rnnum,(float)HithfNS.size());
453  if(HithfNS.size()!= 1728) {
454  cout<<" HF problem "<<rnnum<<" "<<HithfNS.size()<<endl;
455  return;
456  }
457 
458  edm::Handle<HFRecHitCollection> hfMB;
459  iEvent.getByLabel(hfrecoMB, hfMB);
460 
461  if(!hfMB.isValid()){
462  LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
463  cout<<" No HBHE MB"<<endl;
464  return ;
465  }
466 
467  const HFRecHitCollection HithfMB = *(hfMB.product());
468  cout<<" HF MB size of collection "<<HithfMB.size()<<endl;
469  if(HithfMB.size()!= 1728) {
470  cout<<" HF problem "<<rnnum<<" "<<HithfMB.size()<<endl;
471  return;
472  }
473 
474 
475 
476  for(HBHERecHitCollection::const_iterator hbheItr=HithbheNS.begin(); hbheItr!=HithbheNS.end(); hbheItr++)
477  {
478 // Recalibration of energy
479  float icalconst=1.;
480  DetId mydetid = hbheItr->id().rawId();
481  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
482 
483  HBHERecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
484 
485  double energyhit = aHit.energy();
486 
487  DetId id = (*hbheItr).detid();
488  HcalDetId hid=HcalDetId(id);
489 
490  int mysu = ((hid).rawId()>>25)&0x7;
491  if( hid.ieta() > 0 ) {
492  theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
493  theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
494  theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
495  theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
496 
497  tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
498  tmpNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
499 
500 
501  } else {
502  theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
503  theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
504  theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
505  theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
506 
507 
508  tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
509  tmpNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
510 
511  }
512 
513  if(hid.depth() == 1) {
514  if( hid.ieta() > 0 ) {
515  hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit-noise_pl[hid.iphi()][hid.ieta()]);
516  hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
517  } else {
518  hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit-noise_min[hid.iphi()][abs(hid.ieta())]);
519  hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
520  } // eta><0
521  } // depth=1
522  } // HBHE_NS
523 
524 
525 // Signal part for HB HE
526 
527 
528  for(HBHERecHitCollection::const_iterator hbheItr=HithbheMB.begin(); hbheItr!=HithbheMB.end(); hbheItr++)
529  {
530 // Recalibration of energy
531  float icalconst=1.;
532  DetId mydetid = hbheItr->id().rawId();
533  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
534 
535  HBHERecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
536 
537  double energyhit = aHit.energy();
538 
539  DetId id = (*hbheItr).detid();
540  HcalDetId hid=HcalDetId(id);
541 
542  int mysu = ((hid).rawId()>>25)&0x7;
543  if( hid.ieta() > 0 ) {
544  theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
545  theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
546  theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
547  theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
548  float mydiff = energyhit - tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
549 
550 
551  theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+mydiff;
552  theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(mydiff,2);
553  } else {
554  theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
555  theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
556  theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
557  theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
558 
559 
560  float mydiff = energyhit - tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
561  theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+mydiff;
562  theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(mydiff,2);
563  }
564 
565 
566  if(hid.depth() == 1) {
567  if( hid.ieta() > 0 ) {
568  hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit);
569  hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
570  } else {
571  hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit);
572  hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
573  } // eta><0
574  } // depth=1
575  } // HBHE_MB
576 
577 // HF
578 
579  for(HFRecHitCollection::const_iterator hbheItr=HithfNS.begin(); hbheItr!=HithfNS.end(); hbheItr++)
580  {
581 // Recalibration of energy
582  float icalconst=1.;
583  DetId mydetid = hbheItr->id().rawId();
584  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
585 
586  HFRecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
587 
588  double energyhit = aHit.energy();
589 
590  DetId id = (*hbheItr).detid();
591  HcalDetId hid=HcalDetId(id);
592 
593  int mysu = ((hid).rawId()>>25)&0x7;
594  if( hid.ieta() > 0 ) {
595  theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
596  theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
597  theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
598  theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
599 
600  tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
601  tmpNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
602 
603 
604  } else {
605  theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
606  theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
607  theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
608  theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
609 
610 
611  tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
612  tmpNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
613 
614  }
615 
616  if(hid.depth() == 1) {
617  if( hid.ieta() > 0 ) {
618  hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit-noise_pl[hid.iphi()][hid.ieta()]);
619  hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
620  } else {
621  hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit-noise_min[hid.iphi()][abs(hid.ieta())]);
622  hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
623  } // eta><0
624  } // depth=1
625  } // HBHE_NS
626 
627 
628 // Signal part for HB HE
629 
630  for(HFRecHitCollection::const_iterator hbheItr=HithfMB.begin(); hbheItr!=HithfMB.end(); hbheItr++)
631  {
632 // Recalibration of energy
633  float icalconst=1.;
634  DetId mydetid = hbheItr->id().rawId();
635  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
636 
637  HFRecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
638 
639  double energyhit = aHit.energy();
640 
641  DetId id = (*hbheItr).detid();
642  HcalDetId hid=HcalDetId(id);
643 
644  int mysu = ((hid).rawId()>>25)&0x7;
645  if( hid.ieta() > 0 ) {
646  theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
647  theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
648  theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
649  theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
650 
651 
652  theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit-tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
653  theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] =
654  theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow((energyhit-tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]),2);
655  } else {
656  theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
657  theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
658  theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
659  theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
660 
661  theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit-tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
662  theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] =
663  theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow((energyhit-tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]),2);
664  }
665 
666 
667  if(hid.depth() == 1) {
668  if( hid.ieta() > 0 ) {
669  hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit);
670  hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
671  } else {
672  hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit);
673  hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
674  } // eta><0
675  } // depth=1
676  } // HF_MB
677 
678  std::cout<<" Event is finished "<<std::endl;
679 }
680 }
681 //define this as a plug-in
682 //DEFINE_FWK_MODULE(Analyzer_minbias)
683 
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:516
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
hcalLocalRecoNZS_cff.horecoMB
tuple horecoMB
Definition: hcalLocalRecoNZS_cff.py:10
edm::RunBase::run
RunNumber_t run() const
Definition: RunBase.h:44
Analyzer_minbias.h
edm::SortedCollection::const_iterator
std::vector< T >::const_iterator const_iterator
Definition: SortedCollection.h:81
abs
#define abs(x)
Definition: mlp_lapack.h:159
eta
T eta() const
Definition: Basic3DVectorLD.h:168
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:101
bk::beginJob
void beginJob()
Definition: Breakpoints.cc:15
mathSSE::return
return((rh^lh)&mask)
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
nevent
int nevent
Definition: AMPTWrapper.h:74
HcalDetId::depth
int depth() const
get the tower depth
Definition: HcalDetId.h:42
iEvent
int iEvent
Definition: GenABIO.cc:243
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:611
CaloRecHit::energy
float energy() const
Definition: CaloRecHit.h:19
HcalDetId::ieta
int ieta() const
get the cell ieta
Definition: HcalDetId.h:38
edm::Event::run
RunNumber_t run() const
Definition: Event.h:66
hcalLocalRecoNZS_cff.hfrecoMB
tuple hfrecoMB
Definition: hcalLocalRecoNZS_cff.py:7
j
int j
Definition: DBlmapReader.cc:9
ALCARECOHcalCalMinBias_cff.hbherecoNoise
tuple hbherecoNoise
Definition: ALCARECOHcalCalMinBias_cff.py:10
ALCARECOHcalCalMinBias_cff.hfrecoNoise
tuple hfrecoNoise
Definition: ALCARECOHcalCalMinBias_cff.py:12
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
HcalDetId::iphi
int iphi() const
get the cell iphi
Definition: HcalDetId.h:40
DetId
Definition: DetId.h:20
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:62
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
ALCARECOHcalCalMinBias_cff.horecoNoise
tuple horecoNoise
Definition: ALCARECOHcalCalMinBias_cff.py:16
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
hcalLocalRecoNZS_cff.hbherecoMB
tuple hbherecoMB
Definition: hcalLocalRecoNZS_cff.py:4
edm::Run
Definition: Run.h:32
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5