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