CMS 3D CMS Logo

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
27 
28 #include "TFile.h"
29 #include "TH1.h"
30 #include "TH2.h"
31 #include <fstream>
32 #include <sstream>
33 
38 
39 //
40 // constructors and destructor
41 //
42 namespace cms{
44  // get name of output file with histogramms
45  fOutputFileName = iConfig.getUntrackedParameter<std::string>("HistOutFile");
46  // get names of modules, producing object collections
47 
48  tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInputMB"));
49  tok_ho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInputMB"));
50  tok_hf_ = consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfInputMB"));
51  tok_data_ = consumes<FEDRawDataCollection>(edm::InputTag(iConfig.getParameter<std::string>("InputLabel") ));
52 
53  tok_hbheNoise_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInputNoise"));
54  tok_hoNoise_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInputNoise"));
55  tok_hfNoise_ = consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfInputNoise"));
56 
57  // this was hardcodded..
58  tok_gtRec_ = consumes<L1GlobalTriggerReadoutRecord>(edm::InputTag("gtDigisAlCaMB"));
59  tok_hbheNorm_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
60 
61  theRecalib = iConfig.getParameter<bool>("Recalib");
62 
63  //
64  //
65  for(int i=0; i<73; i++) {
66  for(int j=0; j<43; j++) {
67  noise_min[i][j] = 0.;
68  noise_pl[i][j] = 0.;
69  }
70  }
71  //
72  //
73 
74  }
75 
77 
78  // do anything here that needs to be done at desctruction time
79  // (e.g. close files, deallocate resources etc.)
80 
81  }
82 
83  void Analyzer_minbias::beginRun( const edm::Run& r, const edm::EventSetup& iSetup) {
84  nevent_run = 0;
85  }
86  void Analyzer_minbias::endRun( const edm::Run& r, const edm::EventSetup& iSetup) {
87  edm::LogInfo("AnalyzerMB")<<" Runnumber "<<r.run()<<" Nevents "<<nevent_run;
88  }
89 
91 
92  hOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
93 
94  myTree = new TTree("RecJet","RecJet Tree");
95  myTree->Branch("mydet", &mydet, "mydet/I");
96  myTree->Branch("mysubd", &mysubd, "mysubd/I");
97  myTree->Branch("depth", &depth, "depth/I");
98  myTree->Branch("ieta", &ieta, "ieta/I");
99  myTree->Branch("iphi", &iphi, "iphi/I");
100  myTree->Branch("eta", &eta, "eta/F");
101  myTree->Branch("phi", &phi, "phi/F");
102 
103  myTree->Branch("mom0_MB", &mom0_MB, "mom0_MB/F");
104  myTree->Branch("mom1_MB", &mom1_MB, "mom1_MB/F");
105  myTree->Branch("mom2_MB", &mom2_MB, "mom2_MB/F");
106  myTree->Branch("mom4_MB", &mom4_MB, "mom4_MB/F");
107 
108  myTree->Branch("mom0_Noise", &mom0_Noise, "mom0_Noise/F");
109  myTree->Branch("mom1_Noise", &mom1_Noise, "mom1_Noise/F");
110  myTree->Branch("mom2_Noise", &mom2_Noise, "mom2_Noise/F");
111  myTree->Branch("mom4_Noise", &mom4_Noise, "mom4_Noise/F");
112 
113  myTree->Branch("mom0_Diff", &mom0_Diff, "mom0_Diff/F");
114  myTree->Branch("mom1_Diff", &mom1_Diff, "mom1_Diff/F");
115  myTree->Branch("mom2_Diff", &mom2_Diff, "mom2_Diff/F");
116 
117  myTree->Branch("occup", &occup, "occup/F");
118 
119  edm::LogInfo("AnalyzerMB")<<" Before ordering Histos ";
120 
121  char str0[15];
122  char str1[15];
123 
124  char str10[15];
125  char str11[15];
126 
127  int k=0;
128  nevent = 0;
129  // Size of collections
130 
131  hHBHEsize_vs_run = new TH2F("hHBHEsize_vs_run","hHBHEsize_vs_run",500,111500.,112000.,6101,-100.5,6000.5);
132  hHFsize_vs_run = new TH2F("hHFsize_vs_run","hHFsize_vs_run",500,111500.,112000.,6101,-100.5,6000.5);
133 
134  for(int i=1;i<73;i++){
135  for(int j=1;j<43;j++){
136 
137  meannoise_pl[i][j] = 0.;
138  meannoise_min[i][j] = 0.;
139 
140  // for(int l=1;l<5;l++){
141  k = i*1000+j;
142  sprintf(str0,"mpl%d",k);
143  sprintf(str1,"mmin%d",k);
144 
145  sprintf(str10,"vpl%d",k);
146  sprintf(str11,"vmin%d",k);
147  // edm::LogInfo("AnalyzerMB")<<" "<<i<<" "<<j;
148  if( j < 30 ) {
149  // first order moment
150  hCalo1[i][j] = new TH1F(str0, "h0", 320, -10., 10.);
151  hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
152 
153  // second order moment
154  hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 20.);
155  hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 20.);
156  } else {
157  // HF
158  // first order moment
159  // edm::LogInfo("AnalyzerMB")<<" "<<i<<" "<<j<<" "<<k;
160  if(j < 40) {
161  hCalo1[i][j] = new TH1F(str0, "h0", 320, -10., 10.);
162  hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
163  //
164  // second order moment
165  hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 40.);
166  hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 40.);
167  } else {
168  hCalo1[i][j] = new TH1F(str0,"h0" , 320, -10., 10.);
169  hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
170 
171  // second order moment
172  hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 120.);
173  hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 120.);
174 
175  }
176  } // HE/HF boundary
177  // } // l
178  } // j
179  } // i
180 
181 
182  hbheNoiseE = new TH1F("hbheNoiseE","hbheNoiseE", 320, -10., 10.);
183  hfNoiseE = new TH1F("hfNoiseE","hfNoiseE", 320, -10., 10.);
184  hbheSignalE = new TH1F("hbheSignalE","hbheSignalE", 320, -10., 10.);
185  hfSignalE = new TH1F("hfSignalE","hfSignalE", 320, -10., 10.);
186 
187 
188  edm::LogInfo("AnalyzerMB")<<" After ordering Histos ";
189 
190  std::string ccc = "noise_0.dat";
191 
192  myout_hcal = new std::ofstream(ccc.c_str());
193  if(!myout_hcal) edm::LogInfo("AnalyzerMB") << " Output file not open!!! ";
194 
195  //
196  for (int i=0; i<5;i++) {
197  for (int j=0; j<5;j++) {
198  for (int k=0; k<73;k++) {
199  for (int l=0; l<43;l++) {
200  theMBFillDetMapPl0[i][j][k][l] = 0.;
201  theMBFillDetMapPl1[i][j][k][l] = 0.;
202  theMBFillDetMapPl2[i][j][k][l] = 0.;
203  theMBFillDetMapPl4[i][j][k][l] = 0.;
204 
205  theMBFillDetMapMin0[i][j][k][l] = 0.;
206  theMBFillDetMapMin1[i][j][k][l] = 0.;
207  theMBFillDetMapMin2[i][j][k][l] = 0.;
208  theMBFillDetMapMin4[i][j][k][l] = 0.;
209 
210 
211  theNSFillDetMapPl0[i][j][k][l] = 0.;
212  theNSFillDetMapPl1[i][j][k][l] = 0.;
213  theNSFillDetMapPl2[i][j][k][l] = 0.;
214  theNSFillDetMapPl4[i][j][k][l] = 0.;
215 
216  theNSFillDetMapMin0[i][j][k][l] = 0.;
217  theNSFillDetMapMin1[i][j][k][l] = 0.;
218  theNSFillDetMapMin2[i][j][k][l] = 0.;
219  theNSFillDetMapMin4[i][j][k][l] = 0.;
220 
221  theDFFillDetMapPl0[i][j][k][l] = 0.;
222  theDFFillDetMapPl1[i][j][k][l] = 0.;
223  theDFFillDetMapPl2[i][j][k][l] = 0.;
224  theDFFillDetMapMin0[i][j][k][l] = 0.;
225  theDFFillDetMapMin1[i][j][k][l] = 0.;
226  theDFFillDetMapMin2[i][j][k][l] = 0.;
227  }
228  }
229  }
230  }
231 
232  return ;
233  }
234  //
235  // EndJob
236  //
238  int ii=0;
239 
240  for (int i=1; i<5;i++) {
241  for (int j=1; j<5;j++) {
242  for (int k=1; k<73;k++) {
243  for (int l=1; l<43;l++) {
244  if(theMBFillDetMapPl0[i][j][k][l] > 0) {
245  mom0_MB = theMBFillDetMapPl0[i][j][k][l];
246  mom1_MB = theMBFillDetMapPl1[i][j][k][l];
247  mom2_MB = theMBFillDetMapPl2[i][j][k][l];
248  mom4_MB = theMBFillDetMapPl4[i][j][k][l];
253  mom0_Diff = theDFFillDetMapPl0[i][j][k][l];
254  mom1_Diff = theDFFillDetMapPl1[i][j][k][l];
255  mom2_Diff = theDFFillDetMapPl2[i][j][k][l];
256 
257  mysubd = i;
258  depth = j;
259  ieta = l;
260  iphi = k;
261  edm::LogInfo("AnalyzerMB")<<" Result Plus= "<<mysubd<<" "<<ieta<<" "<<iphi<<" mom0 "<<mom0_MB<<" mom1 "<<mom1_MB<<" mom2 "<<mom2_MB;
262  myTree->Fill();
263  ii++;
264  } // Pl > 0
265 
266 
267  if(theMBFillDetMapMin0[i][j][k][l] > 0) {
268  mom0_MB = theMBFillDetMapMin0[i][j][k][l];
269  mom1_MB = theMBFillDetMapMin1[i][j][k][l];
270  mom2_MB = theMBFillDetMapMin2[i][j][k][l];
271  mom4_MB = theMBFillDetMapMin4[i][j][k][l];
279 
280 
281  mysubd = i;
282  depth = j;
283  ieta = -1*l;
284  iphi = k;
285  edm::LogInfo("AnalyzerMB")<<" Result Minus= "<<mysubd<<" "<<ieta<<" "<<iphi<<" mom0 "<<mom0_MB<<" mom1 "<<mom1_MB<<" mom2 "<<mom2_MB;
286  myTree->Fill();
287  ii++;
288 
289  } // Min>0
290  } // ieta
291  } // iphi
292  } // depth
293  } //subd
294 
295 
296 
297  edm::LogInfo("AnalyzerMB")<<" Number of cells "<<ii;
298 
299  hOutputFile->Write();
300 
301  hOutputFile->cd();
302 
303  myTree->Write();
304 
305  hHBHEsize_vs_run->Write() ;
306  hHFsize_vs_run->Write() ;
307 
308  for(int i=1;i<73;i++){
309  for(int j=1;j<43;j++){
310  hCalo1[i][j]->Write();
311  hCalo2[i][j]->Write();
312  hCalo1mom2[i][j]->Write();
313  hCalo2mom2[i][j]->Write();
314  }
315  }
316 
317  hbheNoiseE->Write() ;
318  hfNoiseE->Write() ;
319  hbheSignalE->Write() ;
320  hfSignalE->Write() ;
321 
322 
323  hOutputFile->Close() ;
324 
325  edm::LogInfo("AnalyzerMB")<<" File is closed ";
326 
327  return ;
328  }
329 
330 
331  //
332  // member functions
333  //
334 
335  // ------------ method called to produce the data ------------
337 
338  edm::LogInfo("AnalyzerMB")<<" Start Analyzer_minbias::analyze "<<nevent;
339  nevent++;
340  nevent_run++;
341  using namespace edm;
342 
343  float rnnum = (float)iEvent.run();
344 
345  std::vector<StableProvenance const*> theProvenance;
346  iEvent.getAllStableProvenance(theProvenance);
347 
348  for(auto const& provenance : theProvenance) {
349  edm::LogInfo("AnalyzerMB")<<" Print all process/modulelabel/product names "
350  <<provenance->processName()<<" , "<<provenance->moduleLabel()<<" , "
351  <<provenance->productInstanceName();
352  }
353  /*
354  edm::Handle<FEDRawDataCollection> rawdata;
355  iEvent.getByToken(tok_data_,rawdata);
356 
357  if (!rawdata.isValid()) {
358  edm::LogInfo("AnalyzerMB")<<" No valid collection ";
359  } else {
360  edm::LogInfo("AnalyzerMB")<<" Valid collection ";
361  int calibType = -1 ; int numEmptyFEDs = 0 ;
362  std::vector<int> calibTypeCounter(8,0) ;
363  for (int i=FEDNumbering::MINHCALFEDID;
364  i<=FEDNumbering::MAXHCALFEDID; i++) {
365  const FEDRawData& fedData = rawdata->FEDData(i) ;
366  if ( fedData.size() < 24 ) numEmptyFEDs++ ;
367  if ( fedData.size() < 24 ) continue ;
368  // int value = ((const HcalDCCHeader*)(fedData.data()))->getCalibType() ;
369  // calibTypeCounter.at(value)++ ; // increment the counter for this calib type
370  }
371  edm::LogInfo("AnalyzerMB")<<" NumFed "<<numEmptyFEDs<<" "<<calibType;
372  }
373  */
374  /*
375  std::vector<edm::Handle<FEDRawDataCollection> > rawdata1;
376  iEvent.getManyByType(rawdata1);
377 
378  for(std::vector<edm::Handle<FEDRawDataCollection> >::const_iterator it = rawdata1.begin();it != rawdata1.end(); it++) {
379 
380  edm::LogInfo("AnalyzerMB")<<" Many by Type product name "<< (*it).provenance()->processName()<<
381  " "<<(*it).provenance()->moduleLabel();
382 
383  if((*it).provenance()->moduleLabel() == "hltHcalCalibrationRaw") {
384  int calibType = -1 ; int numEmptyFEDs = 0 ;
385 
386  for (int i=FEDNumbering::MINHCALFEDID;
387  i<=FEDNumbering::MAXHCALFEDID; i++) {
388  const FEDRawData& fedData = (*it)->FEDData(i) ;
389  edm::LogInfo("AnalyzerMB")<<" FED size "<<fedData.size();
390  if ( fedData.size() < 24 ) numEmptyFEDs++ ;
391  if ( fedData.size() < 24 ) continue ;
392  int value = ((const HcalDCCHeader*)(fedData.data()))->getCalibType() ;
393  edm::LogInfo("AnalyzerMB")<<" Value "<<value;
394  }
395  edm::LogInfo("AnalyzerMB")<<" Many by Type NumFed "<<numEmptyFEDs<<" "<<calibType;
396  }
397  }
398 
399 
400  */
401 
402 
403  // Geometry
404  // edm::ESHandle<CaloGeometry> pG;
405  // iSetup.get<CaloGeometryRecord>().get(pG);
406  // ======
407 
408  /*
409  edm::ESHandle<L1GtTriggerMenu> menuRcd;
410  iSetup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
411  const L1GtTriggerMenu* menu = menuRcd.product();
412  const AlgorithmMap& bitMap = menu->gtAlgorithmMap();
413 
414  edm::Handle<L1GlobalTriggerReadoutRecord> gtRecord;
415  iEvent.getByToken(tok_gtRec_, gtRecord);
416 
417  if (!gtRecord.isValid()) {
418 
419  // LogDebug("L1GlobalTriggerRecordProducer")
420  // << "\n\n Error: no L1GlobalTriggerReadoutRecord found with input tag "
421  // << m_l1GtReadoutRecord
422  // << "\n Returning empty L1GlobalTriggerRecord.\n\n";
423  edm::LogInfo("AnalyzerMB")<<" No L1 trigger record ";
424  } else {
425 
426  const DecisionWord dWord = gtRecord->decisionWord();
427 
428  for (CItAlgo itAlgo = bitMap.begin(); itAlgo != bitMap.end(); itAlgo++) {
429  bool decision=menu->gtAlgorithmResult(itAlgo->first,dWord);
430  if(decision == 1) edm::LogInfo("AnalyzerMB")<<" Trigger "<<itAlgo->first<<" "<<decision;
431  }
432 
433  }
434  */
435 
436  const HcalRespCorrs* myRecalib=0;
437  if( theRecalib ) {
438  // Radek:
439  edm::ESHandle <HcalRespCorrs> recalibCorrs;
440  iSetup.get<HcalRespCorrsRcd>().get("recalibrate",recalibCorrs);
441  myRecalib = recalibCorrs.product();
442  // end
443  } // theRecalib
444 
445  // Noise part for HB HE
446 
447  double tmpNSFillDetMapPl1[5][5][73][43];
448  double tmpNSFillDetMapMin1[5][5][73][43];
449 
450  for (int i=0; i<5;i++) {
451  for (int j=0; j<5;j++) {
452  for (int k=0; k<73;k++) {
453  for (int l=0; l<43;l++) {
454  tmpNSFillDetMapPl1[i][j][k][l] = 0.;
455  tmpNSFillDetMapMin1[i][j][k][l] = 0.;
456  }
457  }
458  }
459  }
460 
462  iEvent.getByToken(tok_hbheNorm_, hbheNormal);
463  if(!hbheNormal.isValid()){
464  edm::LogWarning("AnalyzerMB")<<" hbheNormal failed ";
465  } else {
466  edm::LogInfo("AnalyzerMB")<<" The size of the normal collection "<<hbheNormal->size();
467  }
468 
469 
471  iEvent.getByToken(tok_hbheNoise_, hbheNS);
472 
473 
474  if(!hbheNS.isValid()){
475  edm::LogWarning("AnalyzerMB") << "HcalCalibAlgos: Error! can't get hbhe"
476  << " product! No HBHE MS ";
477  return ;
478  }
479 
480 
481  const HBHERecHitCollection HithbheNS = *(hbheNS.product());
482  edm::LogInfo("AnalyzerMB")<<" HBHE NS size of collection "<<HithbheNS.size();
483  hHBHEsize_vs_run->Fill(rnnum,(float)HithbheNS.size());
484 
485  if(HithbheNS.size()!= 5184) {
486  edm::LogWarning("AnalyzerMB")<<" HBHE problem "<<rnnum<<" "<<HithbheNS.size();
487  // return;
488  }
490  iEvent.getByToken(tok_hbhe_, hbheMB);
491 
492  if(!hbheMB.isValid()){
493  edm::LogWarning("AnalyzerMB")<< "HcalCalibAlgos: Error! can't get hbhe"
494  << " product! No HBHE MB";
495  // return ;
496  }
497 
498  const HBHERecHitCollection HithbheMB = *(hbheMB.product());
499  edm::LogInfo("AnalyzerMB")<<" HBHE MB size of collection "<<HithbheMB.size();
500  if(HithbheMB.size()!= 5184) {
501  edm::LogWarning("AnalyzerMB")<<" HBHE problem "<<rnnum<<" "<<HithbheMB.size();
502  // return;
503  }
504 
506  iEvent.getByToken(tok_hfNoise_, hfNS);
507 
508  if(!hfNS.isValid()){
509  edm::LogWarning("AnalyzerMB") << "HcalCalibAlgos: Error! can't get hf"
510  << " product! No HF NS ";
511  // return ;
512  }
513 
514  const HFRecHitCollection HithfNS = *(hfNS.product());
515  edm::LogInfo("AnalyzerMB")<<" HFE NS size of collection "<<HithfNS.size();
516  hHFsize_vs_run->Fill(rnnum,(float)HithfNS.size());
517  if(HithfNS.size()!= 1728) {
518  edm::LogWarning("AnalyzerMB")<<" HF problem "<<rnnum<<" "<<HithfNS.size();
519  // return;
520  }
521 
523  iEvent.getByToken(tok_hf_, hfMB);
524 
525  if(!hfMB.isValid()){
526  edm::LogWarning("AnalyzerMB") << "HcalCalibAlgos: Error! can't get hf"
527  << " product! No HF MB";
528  // return ;
529  }
530 
531  const HFRecHitCollection HithfMB = *(hfMB.product());
532  edm::LogInfo("AnalyzerMB")<<" HF MB size of collection "<<HithfMB.size();
533  if(HithfMB.size()!= 1728) {
534  edm::LogWarning("AnalyzerMB")<<" HF problem "<<rnnum<<" "<<HithfMB.size();
535  // return;
536  }
537 
538 
539 
540  for(HBHERecHitCollection::const_iterator hbheItr=HithbheNS.begin(); hbheItr!=HithbheNS.end(); hbheItr++) {
541  // Recalibration of energy
542  float icalconst=1.;
543  DetId mydetid = hbheItr->id().rawId();
544  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
545 
546  HBHERecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
547 
548  double energyhit = aHit.energy();
549 
550  DetId id = (*hbheItr).detid();
551  HcalDetId hid=HcalDetId(id);
552 
553 
554 
555  int mysu = ((hid).rawId()>>25)&0x7;
556  if( hid.ieta() > 0 ) {
557  theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
558  theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
559  theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
560  theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
561 
562  tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
563 
564 
565  } else {
566  theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
567  theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
568  theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
569  theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
570 
571 
572  tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
573 
574  }
575 
576  if(hid.depth() == 1) {
577  hbheNoiseE->Fill(energyhit);
578 
579  if(energyhit<-2.) edm::LogInfo("AnalyzerMB")<<" Run "<<rnnum<<" ieta,iphi "<<hid.ieta()<<" "<<hid.iphi()<<energyhit;
580 
581  // if( hid.ieta() > 0 ) {
582  // hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit-noise_pl[hid.iphi()][hid.ieta()]);
583  // hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
584  // } else {
585  // hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit-noise_min[hid.iphi()][abs(hid.ieta())]);
586  // hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
587  // } // eta><0
588 
589  } // depth=1
590 
591 
592  } // HBHE_NS
593 
594 
595  // Signal part for HB HE
596 
597 
598  for(HBHERecHitCollection::const_iterator hbheItr=HithbheMB.begin(); hbheItr!=HithbheMB.end(); hbheItr++) {
599  // Recalibration of energy
600  float icalconst=1.;
601  DetId mydetid = hbheItr->id().rawId();
602  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
603 
604  HBHERecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
605 
606  double energyhit = aHit.energy();
607 
608  DetId id = (*hbheItr).detid();
609  HcalDetId hid=HcalDetId(id);
610 
611  int mysu = ((hid).rawId()>>25)&0x7;
612  if( hid.ieta() > 0 ) {
613  theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] += 1.;
614  theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] += energyhit;
615  theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] += pow(energyhit,2);
616  theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] += pow(energyhit,4);
617  float mydiff = energyhit - tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
618 
619 
620  theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+mydiff;
621  theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(mydiff,2);
622  } else {
623  theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
624  theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
625  theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
626  theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
627 
628 
629  float mydiff = energyhit - tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
630  theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+mydiff;
631  theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(mydiff,2);
632  }
633 
634 
635  if(hid.depth() == 1) {
636 
637  hbheSignalE->Fill(energyhit);
638 
639  if( hid.ieta() > 0 ) {
640  hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit);
641  hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
642  } else {
643  hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit);
644  hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
645  } // eta><0
646 
647  } // depth=1
648 
649 
650  } // HBHE_MB
651 
652  // HF
653 
654  for(HFRecHitCollection::const_iterator hbheItr=HithfNS.begin(); hbheItr!=HithfNS.end(); hbheItr++) {
655  // Recalibration of energy
656  float icalconst=1.;
657  DetId mydetid = hbheItr->id().rawId();
658  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
659 
660  HFRecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
661 
662  double energyhit = aHit.energy();
663  //
664  // Remove PMT hits
665  //
666  DetId id = (*hbheItr).detid();
667  HcalDetId hid=HcalDetId(id);
668  // theGeometry->getGeometry(detId)->getPosition()).eta()
669 
670  // edm::LogInfo("AnalyzerMB")<<hid.ieta()<<" "<<hid.iphi()<<" "<<hid.depth()<<" "<<(pG->getGeometry(id)->getPosition()).eta()<<" "<<(pG->getGeometry(id)->getPosition()).phi()
671  // <<" "<<(pG->getGeometry(id)->getPosition()).perp()<<" "<<(pG->getGeometry(id)->getPosition()).z();
672 
673  if(fabs(energyhit) > 40. ) continue;
674 
675  int mysu = ((hid).rawId()>>25)&0x7;
676  if( hid.ieta() > 0 ) {
677  theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
678  theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
679  theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
680  theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
681 
682  tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
683 
684 
685  } else {
686  theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
687  theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
688  theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
689  theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
690 
691 
692  tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
693 
694  }
695 
696  if(hid.depth() == 1) {
697  hfNoiseE->Fill(energyhit);
698 
699  //if( hid.ieta() > 0 ) {
700  // hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit-noise_pl[hid.iphi()][hid.ieta()]);
701  // hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
702  //} else {
703  // hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit-noise_min[hid.iphi()][abs(hid.ieta())]);
704  // hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
705  //} // eta><0
706 
707  } // depth=1
708 
709  } // HBHE_NS
710 
711 
712  // Signal part for HB HE
713 
714  for(HFRecHitCollection::const_iterator hbheItr=HithfMB.begin(); hbheItr!=HithfMB.end(); hbheItr++) {
715  // Recalibration of energy
716  float icalconst=1.;
717  DetId mydetid = hbheItr->id().rawId();
718  if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
719 
720  HFRecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
721 
722  double energyhit = aHit.energy();
723  //
724  // Remove PMT hits
725  //
726  if(fabs(energyhit) > 40. ) continue;
727 
728  DetId id = (*hbheItr).detid();
729  HcalDetId hid=HcalDetId(id);
730 
731  int mysu = ((hid).rawId()>>25)&0x7;
732  if( hid.ieta() > 0 ) {
733  theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
734  theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
735  theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
736  theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
737 
738 
739  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()];
740  theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] =
741  theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow((energyhit-tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]),2);
742  } else {
743  theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
744  theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
745  theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
746  theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
747 
748  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()];
749  theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] =
750  theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow((energyhit-tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]),2);
751  }
752 
753 
754  if(hid.depth() == 1) {
755  hfSignalE->Fill(energyhit);
756 
757  if( hid.ieta() > 0 ) {
758  hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit);
759  hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
760  } else {
761  hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit);
762  hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
763  } // eta><0
764 
765  } // depth=1
766 
767  } // HF_MB
768 
769  edm::LogInfo("AnalyzerMB")<<" Event is finished ";
770  }
771 }
772 
774 
776 
double theMBFillDetMapMin0[5][5][73][43]
double theMBFillDetMapMin2[5][5][73][43]
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
double theDFFillDetMapPl2[5][5][73][43]
double theMBFillDetMapMin1[5][5][73][43]
RunNumber_t run() const
Definition: RunBase.h:40
double theNSFillDetMapPl4[5][5][73][43]
double theNSFillDetMapPl0[5][5][73][43]
double theNSFillDetMapMin0[5][5][73][43]
double theMBFillDetMapPl0[5][5][73][43]
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:460
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
TH1F * hCalo2mom2[73][43]
edm::EDGetTokenT< HORecHitCollection > tok_hoNoise_
std::vector< HBHERecHit >::const_iterator const_iterator
double theDFFillDetMapPl1[5][5][73][43]
virtual void beginRun(const edm::Run &r, const edm::EventSetup &iSetup)
std::ofstream * myout_hcal
edm::EDGetTokenT< HFRecHitCollection > tok_hf_
virtual void endRun(const edm::Run &r, const edm::EventSetup &iSetup)
edm::EDGetTokenT< L1GlobalTriggerReadoutRecord > tok_gtRec_
const Item * getValues(DetId fId, bool throwOnFail=true) const
double meannoise_pl[73][43]
edm::EDGetTokenT< HBHERecHitCollection > tok_hbheNoise_
double theNSFillDetMapPl2[5][5][73][43]
double theDFFillDetMapMin2[5][5][73][43]
double noise_min[73][43]
TH1F * hCalo1mom2[73][43]
edm::EDGetTokenT< HFRecHitCollection > tok_hfNoise_
double theNSFillDetMapPl1[5][5][73][43]
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
double theMBFillDetMapPl1[5][5][73][43]
int depth() const
get the tower depth
Definition: HcalDetId.cc:108
Analyzer_minbias(const edm::ParameterSet &)
int iEvent
Definition: GenABIO.cc:230
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
float energy() const
Definition: CaloRecHit.h:17
double theDFFillDetMapPl0[5][5][73][43]
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
int ieta() const
get the cell ieta
Definition: HcalDetId.h:56
RunNumber_t run() const
Definition: Event.h:94
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double theMBFillDetMapPl4[5][5][73][43]
double theNSFillDetMapMin1[5][5][73][43]
double theNSFillDetMapMin4[5][5][73][43]
bool isValid() const
Definition: HandleBase.h:74
double theDFFillDetMapMin0[5][5][73][43]
double theMBFillDetMapPl2[5][5][73][43]
edm::EDGetTokenT< HBHERecHitCollection > tok_hbheNorm_
ii
Definition: cuy.py:588
int k[5][pyjets_maxn]
const_iterator end() const
int iphi() const
get the cell iphi
Definition: HcalDetId.cc:103
Definition: DetId.h:18
virtual void analyze(const edm::Event &, const edm::EventSetup &)
T const * product() const
Definition: Handle.h:81
edm::EDGetTokenT< FEDRawDataCollection > tok_data_
const T & get() const
Definition: EventSetup.h:56
double theMBFillDetMapMin4[5][5][73][43]
return(e1-e2)*(e1-e2)+dp *dp
double meannoise_min[73][43]
HLT enums.
float getValue() const
Definition: HcalRespCorr.h:20
size_type size() const
edm::EDGetTokenT< HORecHitCollection > tok_ho_
double theNSFillDetMapMin2[5][5][73][43]
double theDFFillDetMapMin1[5][5][73][43]
T const * product() const
Definition: ESHandle.h:86
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
const_iterator begin() const
Definition: Run.h:42
void getAllStableProvenance(std::vector< StableProvenance const * > &provenances) const
Definition: Event.cc:96