Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:59 

0001 // system include files
0002 #include <memory>
0003 #include <string>
0004 #include <iostream>
0005 #include <fstream>
0006 #include <sstream>
0007 #include <vector>
0008 #include <map>
0009 
0010 // user include files
0011 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0012 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
0013 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0014 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
0015 
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0018 
0019 #include "DataFormats/Common/interface/Ref.h"
0020 #include "DataFormats/Common/interface/TriggerResults.h"
0021 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0022 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0023 #include "DataFormats/Candidate/interface/Candidate.h"
0024 #include "DataFormats/DetId/interface/DetId.h"
0025 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0026 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0027 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0028 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0029 #include "DataFormats/HcalDigi/interface/HcalCalibrationEventTypes.h"
0030 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0031 #include "DataFormats/JetReco/interface/Jet.h"
0032 #include "DataFormats/JetReco/interface/CaloJet.h"
0033 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0034 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
0035 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0036 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0037 #include "DataFormats/L1GlobalTrigger/interface/L1GtfeWord.h"
0038 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0039 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerRecord.h"
0040 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0041 #include "DataFormats/Provenance/interface/Provenance.h"
0042 
0043 #include "EventFilter/HcalRawToDigi/interface/HcalDCCHeader.h"
0044 
0045 #include "FWCore/Common/interface/TriggerNames.h"
0046 #include "FWCore/Framework/interface/Frameworkfwd.h"
0047 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0048 #include "FWCore/Framework/interface/Event.h"
0049 #include "FWCore/Framework/interface/EventSetup.h"
0050 #include "FWCore/Framework/interface/Run.h"
0051 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0052 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0053 
0054 #include "TFile.h"
0055 #include "TH1.h"
0056 #include "TH2.h"
0057 #include "TTree.h"
0058 
0059 namespace HcalMinbias {}
0060 
0061 // constructors and destructor
0062 class AnalyzerMinbias : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0063 public:
0064   explicit AnalyzerMinbias(const edm::ParameterSet&);
0065   ~AnalyzerMinbias() override = default;
0066 
0067   void analyze(const edm::Event&, const edm::EventSetup&) override;
0068   void beginJob() override;
0069   void endJob() override;
0070 
0071 private:
0072   void analyzeHcal(const HcalRespCorrs* myRecalib,
0073                    const HBHERecHitCollection& HithbheNS,
0074                    const HBHERecHitCollection& HithbheMB,
0075                    const HFRecHitCollection& HithfNS,
0076                    const HFRecHitCollection& HithfMB,
0077                    int algoBit,
0078                    bool fill);
0079 
0080   struct myInfo {
0081     double theMB0, theMB1, theMB2, theMB3, theMB4;
0082     double theNS0, theNS1, theNS2, theNS3, theNS4;
0083     double theDif0, theDif1, theDif2, runcheck;
0084     myInfo() {
0085       theMB0 = theMB1 = theMB2 = theMB3 = theMB4 = 0;
0086       theNS0 = theNS1 = theNS2 = theNS3 = theNS4 = 0;
0087       theDif0 = theDif1 = theDif2 = runcheck = 0;
0088     }
0089   };
0090 
0091   // ----------member data ---------------------------
0092   const std::string fOutputFileName;
0093   const bool runNZS_, theRecalib_, ignoreL1_;
0094   std::string hcalfile_;
0095   std::ofstream* myout_hcal;
0096   TFile* hOutputFile_;
0097   TTree* myTree_;
0098   TH1D *h_Noise[4], *h_Signal[4];
0099   double rnnum_;
0100 
0101   // Root tree members
0102   double rnnumber;
0103   int mydet, mysubd, depth, iphi, ieta, cells, trigbit;
0104   float phi, eta;
0105   float mom0_MB, mom1_MB, mom2_MB, mom3_MB, mom4_MB, occup;
0106   float mom0_Noise, mom1_Noise, mom2_Noise, mom3_Noise, mom4_Noise;
0107   float mom0_Diff, mom1_Diff, mom2_Diff, mom3_Diff, mom4_Diff;
0108 
0109   std::map<std::pair<int, HcalDetId>, myInfo> myMap_;
0110   const edm::EDGetTokenT<HBHERecHitCollection> tok_hbherecoMB_, tok_hbherecoNoise_;
0111   const edm::EDGetTokenT<HORecHitCollection> tok_horecoMB_, tok_horecoNoise_;
0112   const edm::EDGetTokenT<HFRecHitCollection> tok_hfrecoMB_, tok_hfrecoNoise_;
0113   const edm::EDGetTokenT<HBHERecHitCollection> tok_hbheNormal_;
0114   const edm::EDGetTokenT<L1GlobalTriggerObjectMapRecord> tok_hltL1GtMap_;
0115   const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_respCorr_;
0116 };
0117 
0118 AnalyzerMinbias::AnalyzerMinbias(const edm::ParameterSet& iConfig)
0119     : fOutputFileName(iConfig.getUntrackedParameter<std::string>("HistOutFile")),
0120       runNZS_(iConfig.getUntrackedParameter<bool>("RunNZS", true)),
0121       theRecalib_(iConfig.getParameter<bool>("Recalib")),
0122       ignoreL1_(iConfig.getUntrackedParameter<bool>("IgnoreL1", true)),
0123       tok_hbherecoMB_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInputMB"))),
0124       tok_hbherecoNoise_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInputNoise"))),
0125       tok_horecoMB_(consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInputMB"))),
0126       tok_horecoNoise_(consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInputNoise"))),
0127       tok_hfrecoMB_(consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfInputMB"))),
0128       tok_hfrecoNoise_(consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfInputNoise"))),
0129       tok_hbheNormal_(consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"))),
0130       tok_hltL1GtMap_(consumes<L1GlobalTriggerObjectMapRecord>(edm::InputTag("hltL1GtObjectMap"))),
0131       tok_respCorr_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()) {
0132   usesResource(TFileService::kSharedResource);
0133   // get name of output file with histogramms
0134   // get token names of modules, producing object collections
0135 }
0136 
0137 void AnalyzerMinbias::beginJob() {
0138   std::string det[4] = {"HB", "HE", "HO", "HF"};
0139   char name[80], title[80];
0140   for (int subd = 0; subd < 4; ++subd) {
0141     sprintf(name, "Noise_%s", det[subd].c_str());
0142     sprintf(title, "Energy Distribution for Noise in %s", det[subd].c_str());
0143     h_Noise[subd] = new TH1D(name, title, 100, -10., 10.);
0144     sprintf(name, "Signal_%s", det[subd].c_str());
0145     sprintf(title, "Energy Distribution for Signal in %s", det[subd].c_str());
0146     h_Signal[subd] = new TH1D(name, title, 100, -10., 10.);
0147   }
0148 
0149   edm::Service<TFileService> fs;
0150   hOutputFile_ = new TFile(fOutputFileName.c_str(), "RECREATE");
0151   myTree_ = fs->make<TTree>("RecJet", "RecJet Tree");
0152   myTree_->Branch("mydet", &mydet, "mydet/I");
0153   myTree_->Branch("mysubd", &mysubd, "mysubd/I");
0154   myTree_->Branch("cells", &cells, "cells");
0155   myTree_->Branch("depth", &depth, "depth/I");
0156   myTree_->Branch("ieta", &ieta, "ieta/I");
0157   myTree_->Branch("iphi", &iphi, "iphi/I");
0158   myTree_->Branch("eta", &eta, "eta/F");
0159   myTree_->Branch("phi", &phi, "phi/F");
0160   myTree_->Branch("mom0_MB", &mom0_MB, "mom0_MB/F");
0161   myTree_->Branch("mom1_MB", &mom1_MB, "mom1_MB/F");
0162   myTree_->Branch("mom2_MB", &mom2_MB, "mom2_MB/F");
0163   myTree_->Branch("mom3_MB", &mom3_MB, "mom3_MB/F");
0164   myTree_->Branch("mom4_MB", &mom4_MB, "mom4_MB/F");
0165   myTree_->Branch("mom0_Noise", &mom0_Noise, "mom0_Noise/F");
0166   myTree_->Branch("mom1_Noise", &mom1_Noise, "mom1_Noise/F");
0167   myTree_->Branch("mom2_Noise", &mom2_Noise, "mom2_Noise/F");
0168   myTree_->Branch("mom3_Noise", &mom3_Noise, "mom3_Noise/F");
0169   myTree_->Branch("mom4_Noise", &mom4_Noise, "mom4_Noise/F");
0170   myTree_->Branch("mom0_Diff", &mom0_Diff, "mom0_Diff/F");
0171   myTree_->Branch("mom1_Diff", &mom1_Diff, "mom1_Diff/F");
0172   myTree_->Branch("mom2_Diff", &mom2_Diff, "mom2_Diff/F");
0173   myTree_->Branch("occup", &occup, "occup/F");
0174   myTree_->Branch("trigbit", &trigbit, "trigbit/I");
0175   myTree_->Branch("rnnumber", &rnnumber, "rnnumber/D");
0176 
0177   myMap_.clear();
0178 }
0179 
0180 //  EndJob
0181 //
0182 void AnalyzerMinbias::endJob() {
0183   int ii = 0;
0184   for (std::map<std::pair<int, HcalDetId>, myInfo>::const_iterator itr = myMap_.begin(); itr != myMap_.end(); ++itr) {
0185 #ifdef EDM_ML_DEBUG
0186     edm::LogVerbatim("AnalyzerMinimumBias") << "Fired trigger bit number " << itr->first.first;
0187 #endif
0188     myInfo info = itr->second;
0189     if (info.theMB0 > 0) {
0190       mom0_MB = info.theMB0;
0191       mom1_MB = info.theMB1;
0192       mom2_MB = info.theMB2;
0193       mom3_MB = info.theMB3;
0194       mom4_MB = info.theMB4;
0195       mom0_Noise = info.theNS0;
0196       mom1_Noise = info.theNS1;
0197       mom2_Noise = info.theNS2;
0198       mom3_Noise = info.theNS3;
0199       mom4_Noise = info.theNS4;
0200       mom0_Diff = info.theDif0;
0201       mom1_Diff = info.theDif1;
0202       mom2_Diff = info.theDif2;
0203       rnnumber = info.runcheck;
0204       trigbit = itr->first.first;
0205       mysubd = itr->first.second.subdet();
0206       depth = itr->first.second.depth();
0207       ieta = itr->first.second.ieta();
0208       iphi = itr->first.second.iphi();
0209 #ifdef EDM_ML_DEBUG
0210       edm::LogVerbatim("AnalyzerMinimumBias")
0211           << " Result=  " << trigbit << " " << mysubd << " " << ieta << " " << iphi << " mom0  " << mom0_MB << " mom1 "
0212           << mom1_MB << " mom2 " << mom2_MB << " mom3 " << mom3_MB << " mom4 " << mom4_MB << " mom0_Noise "
0213           << mom0_Noise << " mom1_Noise " << mom1_Noise << " mom2_Noise " << mom2_Noise << " mom3_Noise " << mom3_Noise
0214           << " mom4_Noise " << mom4_Noise << " mom0_Diff " << mom0_Diff << " mom1_Diff " << mom1_Diff << " mom2_Diff "
0215           << mom2_Diff;
0216 #endif
0217       myTree_->Fill();
0218       ii++;
0219     }
0220   }
0221   cells = ii;
0222 #ifdef EDM_ML_DEBUG
0223   edm::LogVerbatim("AnalyzerMinimumBias") << "cells " << cells;
0224 #endif
0225   hOutputFile_->Write();
0226   hOutputFile_->cd();
0227   myTree_->Write();
0228   for (int i = 0; i < 4; i++) {
0229     h_Noise[i]->Write();
0230     h_Signal[i]->Write();
0231   }
0232   hOutputFile_->Close();
0233 }
0234 
0235 //
0236 // member functions
0237 //
0238 
0239 // ------------ method called to produce the data  ------------
0240 
0241 void AnalyzerMinbias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0242   rnnum_ = (float)iEvent.run();
0243   const HcalRespCorrs* myRecalib = nullptr;
0244   if (theRecalib_) {
0245     myRecalib = &iSetup.getData(tok_respCorr_);
0246   }  // theRecalib
0247 
0248   const edm::Handle<HBHERecHitCollection> hbheNormal = iEvent.getHandle(tok_hbheNormal_);
0249   if (!hbheNormal.isValid()) {
0250     edm::LogVerbatim("AnalyzerMinimumBias") << " hbheNormal failed";
0251   } else {
0252     edm::LogVerbatim("AnalyzerMinimumBias") << " The size of the normal collection " << hbheNormal->size();
0253   }
0254 
0255   const edm::Handle<HBHERecHitCollection> hbheNS = iEvent.getHandle(tok_hbherecoNoise_);
0256   if (!hbheNS.isValid()) {
0257     edm::LogWarning("AnalyzerMinimumBias") << "HcalCalibAlgos: Error! can't get hbheNoise product!";
0258     return;
0259   }
0260   const HBHERecHitCollection HithbheNS = *(hbheNS.product());
0261   edm::LogVerbatim("AnalyzerMinimumBias") << "HBHE NS size of collection " << HithbheNS.size();
0262   if (runNZS_ && HithbheNS.size() != 5184) {
0263     edm::LogWarning("AnalyzerMinimumBias") << "HBHE NS problem " << rnnum_ << " size " << HithbheNS.size();
0264     return;
0265   }
0266 
0267   const edm::Handle<HBHERecHitCollection> hbheMB = iEvent.getHandle(tok_hbherecoMB_);
0268   if (!hbheMB.isValid()) {
0269     edm::LogWarning("AnalyzerMinimumBias") << "HcalCalibAlgos: Error! can't get hbhe product!";
0270     return;
0271   }
0272   const HBHERecHitCollection HithbheMB = *(hbheMB.product());
0273   edm::LogVerbatim("AnalyzerMinimumBias") << "HBHE MB size of collection " << HithbheMB.size();
0274   if (runNZS_ && HithbheMB.size() != 5184) {
0275     edm::LogWarning("AnalyzerMinimumBias") << "HBHE problem " << rnnum_ << " size " << HithbheMB.size();
0276     return;
0277   }
0278 
0279   const edm::Handle<HFRecHitCollection> hfNS = iEvent.getHandle(tok_hfrecoNoise_);
0280   if (!hfNS.isValid()) {
0281     edm::LogWarning("AnalyzerMinimumBias") << "HcalCalibAlgos: Error! can't get hfNoise product!";
0282     return;
0283   }
0284   const HFRecHitCollection HithfNS = *(hfNS.product());
0285   edm::LogVerbatim("AnalyzerMinimumBias") << "HF NS size of collection " << HithfNS.size();
0286   if (runNZS_ && HithfNS.size() != 1728) {
0287     edm::LogWarning("AnalyzerMinimumBias") << "HF NS problem " << rnnum_ << " size " << HithfNS.size();
0288     return;
0289   }
0290 
0291   const edm::Handle<HFRecHitCollection> hfMB = iEvent.getHandle(tok_hfrecoMB_);
0292   if (!hfMB.isValid()) {
0293     edm::LogWarning("AnalyzerMinimumBias") << "HcalCalibAlgos: Error! can't get hf product!";
0294     return;
0295   }
0296   const HFRecHitCollection HithfMB = *(hfMB.product());
0297   edm::LogVerbatim("AnalyzerMinimumBias") << "HF MB size of collection " << HithfMB.size();
0298   if (runNZS_ && HithfMB.size() != 1728) {
0299     edm::LogWarning("AnalyzerMinimumBias") << "HF problem " << rnnum_ << " size " << HithfMB.size();
0300     return;
0301   }
0302 
0303   if (ignoreL1_) {
0304     analyzeHcal(myRecalib, HithbheNS, HithbheMB, HithfNS, HithfMB, 1, true);
0305   } else {
0306     const edm::Handle<L1GlobalTriggerObjectMapRecord> gtObjectMapRecord = iEvent.getHandle(tok_hltL1GtMap_);
0307     if (gtObjectMapRecord.isValid()) {
0308       const std::vector<L1GlobalTriggerObjectMap>& objMapVec = gtObjectMapRecord->gtObjectMap();
0309       int ii(0);
0310       bool ok(false), fill(true);
0311       for (std::vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin(); itMap != objMapVec.end();
0312            ++itMap, ++ii) {
0313         bool resultGt = (*itMap).algoGtlResult();
0314         if (resultGt == 1) {
0315           ok = true;
0316           int algoBit = (*itMap).algoBitNumber();
0317           analyzeHcal(myRecalib, HithbheNS, HithbheMB, HithfNS, HithfMB, algoBit, fill);
0318           fill = false;
0319           std::string algoNameStr = (*itMap).algoName();
0320 #ifdef EDM_ML_DEBUG
0321           edm::LogVerbatim("AnalyzerMinimumBias")
0322               << "Trigger[" << ii << "] " << algoNameStr << " bit " << algoBit << " entered";
0323 #endif
0324         }
0325       }
0326       if (!ok)
0327         edm::LogVerbatim("AnalyzerMinimumBias") << "No passed L1 Triggers";
0328     }
0329   }
0330 }
0331 
0332 void AnalyzerMinbias::analyzeHcal(const HcalRespCorrs* myRecalib,
0333                                   const HBHERecHitCollection& HithbheNS,
0334                                   const HBHERecHitCollection& HithbheMB,
0335                                   const HFRecHitCollection& HithfNS,
0336                                   const HFRecHitCollection& HithfMB,
0337                                   int algoBit,
0338                                   bool fill) {
0339   // Noise part for HB HE
0340   std::map<std::pair<int, HcalDetId>, myInfo> tmpMap;
0341   tmpMap.clear();
0342 
0343   for (HBHERecHitCollection::const_iterator hbheItr = HithbheNS.begin(); hbheItr != HithbheNS.end(); hbheItr++) {
0344     // Recalibration of energy
0345     float icalconst = 1.;
0346     DetId mydetid = hbheItr->id().rawId();
0347     if (theRecalib_)
0348       icalconst = myRecalib->getValues(mydetid)->getValue();
0349 
0350     HBHERecHit aHit(hbheItr->id(), hbheItr->energy() * icalconst, hbheItr->time());
0351     double energyhit = aHit.energy();
0352 
0353     DetId id = (*hbheItr).detid();
0354     HcalDetId hid = HcalDetId(id);
0355     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0356     if (itr1 == myMap_.end()) {
0357       myInfo info;
0358       myMap_[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0359       itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0360     }
0361     itr1->second.theNS0++;
0362     itr1->second.theNS1 += energyhit;
0363     itr1->second.theNS2 += (energyhit * energyhit);
0364     itr1->second.theNS3 += (energyhit * energyhit * energyhit);
0365     itr1->second.theNS4 += (energyhit * energyhit * energyhit * energyhit);
0366     itr1->second.runcheck = rnnum_;
0367     if (fill)
0368       h_Noise[hid.subdet() - 1]->Fill(energyhit);
0369 
0370     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr2 = tmpMap.find(std::pair<int, HcalDetId>(algoBit, hid));
0371     if (itr2 == tmpMap.end()) {
0372       myInfo info;
0373       tmpMap[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0374       itr2 = tmpMap.find(std::pair<int, HcalDetId>(algoBit, hid));
0375     }
0376     itr2->second.theNS0++;
0377     itr2->second.theNS1 += energyhit;
0378     itr2->second.theNS2 += (energyhit * energyhit);
0379     itr2->second.theNS3 += (energyhit * energyhit * energyhit);
0380     itr2->second.theNS4 += (energyhit * energyhit * energyhit * energyhit);
0381     itr2->second.runcheck = rnnum_;
0382 
0383   }  // HBHE_NS
0384 
0385   // Signal part for HB HE
0386 
0387   for (HBHERecHitCollection::const_iterator hbheItr = HithbheMB.begin(); hbheItr != HithbheMB.end(); hbheItr++) {
0388     // Recalibration of energy
0389     float icalconst = 1.;
0390     DetId mydetid = hbheItr->id().rawId();
0391     if (theRecalib_)
0392       icalconst = myRecalib->getValues(mydetid)->getValue();
0393 
0394     HBHERecHit aHit(hbheItr->id(), hbheItr->energy() * icalconst, hbheItr->time());
0395     double energyhit = aHit.energy();
0396 
0397     DetId id = (*hbheItr).detid();
0398     HcalDetId hid = HcalDetId(id);
0399 
0400     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0401     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr2 = tmpMap.find(std::pair<int, HcalDetId>(algoBit, hid));
0402 
0403     if (itr1 == myMap_.end()) {
0404       myInfo info;
0405       myMap_[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0406       itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0407     }
0408     itr1->second.theMB0++;
0409     itr1->second.theDif0 = 0;
0410     itr1->second.theMB1 += energyhit;
0411     itr1->second.theMB2 += (energyhit * energyhit);
0412     itr1->second.theMB3 += (energyhit * energyhit * energyhit);
0413     itr1->second.theMB4 += (energyhit * energyhit * energyhit * energyhit);
0414     itr1->second.runcheck = rnnum_;
0415     float mydiff = 0.0;
0416     if (itr2 != tmpMap.end()) {
0417       mydiff = energyhit - (itr2->second.theNS1);
0418       itr1->second.theDif0++;
0419       itr1->second.theDif1 += mydiff;
0420       itr1->second.theDif2 += (mydiff * mydiff);
0421       if (fill)
0422         h_Signal[hid.subdet() - 1]->Fill(mydiff);
0423     }
0424   }  // HBHE_MB
0425 
0426   // HF
0427 
0428   for (HFRecHitCollection::const_iterator hbheItr = HithfNS.begin(); hbheItr != HithfNS.end(); hbheItr++) {
0429     // Recalibration of energy
0430     float icalconst = 1.;
0431     DetId mydetid = hbheItr->id().rawId();
0432     if (theRecalib_)
0433       icalconst = myRecalib->getValues(mydetid)->getValue();
0434 
0435     HFRecHit aHit(hbheItr->id(), hbheItr->energy() * icalconst, hbheItr->time());
0436     double energyhit = aHit.energy();
0437     // Remove PMT hits
0438     if (fabs(energyhit) > 40.)
0439       continue;
0440     DetId id = (*hbheItr).detid();
0441     HcalDetId hid = HcalDetId(id);
0442 
0443     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0444 
0445     if (itr1 == myMap_.end()) {
0446       myInfo info;
0447       myMap_[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0448       itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0449     }
0450     itr1->second.theNS0++;
0451     itr1->second.theNS1 += energyhit;
0452     itr1->second.theNS2 += (energyhit * energyhit);
0453     itr1->second.theNS3 += (energyhit * energyhit * energyhit);
0454     itr1->second.theNS4 += (energyhit * energyhit * energyhit * energyhit);
0455     itr1->second.runcheck = rnnum_;
0456     if (fill)
0457       h_Noise[hid.subdet() - 1]->Fill(energyhit);
0458 
0459     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr2 = tmpMap.find(std::pair<int, HcalDetId>(algoBit, hid));
0460     if (itr2 == tmpMap.end()) {
0461       myInfo info;
0462       tmpMap[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0463       itr2 = tmpMap.find(std::pair<int, HcalDetId>(algoBit, hid));
0464     }
0465     itr2->second.theNS0++;
0466     itr2->second.theNS1 += energyhit;
0467     itr2->second.theNS2 += (energyhit * energyhit);
0468     itr2->second.theNS3 += (energyhit * energyhit * energyhit);
0469     itr2->second.theNS4 += (energyhit * energyhit * energyhit * energyhit);
0470     itr2->second.runcheck = rnnum_;
0471 
0472   }  // HF_NS
0473 
0474   // Signal part for HF
0475 
0476   for (HFRecHitCollection::const_iterator hbheItr = HithfMB.begin(); hbheItr != HithfMB.end(); hbheItr++) {
0477     // Recalibration of energy
0478     float icalconst = 1.;
0479     DetId mydetid = hbheItr->id().rawId();
0480     if (theRecalib_)
0481       icalconst = myRecalib->getValues(mydetid)->getValue();
0482     HFRecHit aHit(hbheItr->id(), hbheItr->energy() * icalconst, hbheItr->time());
0483 
0484     double energyhit = aHit.energy();
0485     // Remove PMT hits
0486     if (fabs(energyhit) > 40.)
0487       continue;
0488 
0489     DetId id = (*hbheItr).detid();
0490     HcalDetId hid = HcalDetId(id);
0491 
0492     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0493     std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr2 = tmpMap.find(std::pair<int, HcalDetId>(algoBit, hid));
0494 
0495     if (itr1 == myMap_.end()) {
0496       myInfo info;
0497       myMap_[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0498       itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0499     }
0500     itr1->second.theMB0++;
0501     itr1->second.theDif0 = 0;
0502     itr1->second.theMB1 += energyhit;
0503     itr1->second.theMB2 += (energyhit * energyhit);
0504     itr1->second.theMB3 += (energyhit * energyhit * energyhit);
0505     itr1->second.theMB4 += (energyhit * energyhit * energyhit * energyhit);
0506     itr1->second.runcheck = rnnum_;
0507     float mydiff = 0.0;
0508     if (itr2 != tmpMap.end()) {
0509       mydiff = energyhit - (itr2->second.theNS1);
0510       itr1->second.theDif0++;
0511       itr1->second.theDif1 += mydiff;
0512       itr1->second.theDif2 += (mydiff * mydiff);
0513       if (fill)
0514         h_Signal[hid.subdet() - 1]->Fill(mydiff);
0515     }
0516   }
0517 }
0518 
0519 //define this as a plug-in
0520 #include "FWCore/Framework/interface/MakerMacros.h"
0521 DEFINE_FWK_MODULE(AnalyzerMinbias);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!