File indexing completed on 2024-04-06 12:29:29
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #include <memory>
0019 #include <utility>
0020
0021
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "DataFormats/EcalDigi/interface/EcalTriggerPrimitiveDigi.h"
0024 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0025 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0026 #include "CalibCalorimetry/EcalTPGTools/interface/EcalTPGScale.h"
0027 #include "EcalTrigPrimAnalyzer.h"
0028
0029 #include <TMath.h>
0030
0031 class CaloSubdetectorGeometry;
0032
0033 EcalTrigPrimAnalyzer::EcalTrigPrimAnalyzer(const edm::ParameterSet &iConfig)
0034 : recHits_(iConfig.getParameter<bool>("AnalyzeRecHits")),
0035 label_(iConfig.getParameter<edm::InputTag>("inputTP")),
0036 rechits_labelEB_(iConfig.getParameter<edm::InputTag>("inputRecHitsEB")),
0037 rechits_labelEE_(iConfig.getParameter<edm::InputTag>("inputRecHitsEE")),
0038 tpToken_(consumes<EcalTrigPrimDigiCollection>(label_)),
0039 ebToken_(consumes<EcalRecHitCollection>(rechits_labelEB_)),
0040 eeToken_(consumes<EcalRecHitCollection>(rechits_labelEE_)),
0041 tokens_(consumesCollector()) {
0042 usesResource(TFileService::kSharedResource);
0043
0044 ecal_parts_.push_back("Barrel");
0045 ecal_parts_.push_back("Endcap");
0046
0047 edm::Service<TFileService> fs;
0048 tree_ = fs->make<TTree>("TPGtree", "TPGtree");
0049 tree_->Branch("iphi", &iphi_, "iphi/I");
0050 tree_->Branch("ieta", &ieta_, "ieta/I");
0051 tree_->Branch("eRec", &eRec_, "eRec/F");
0052 tree_->Branch("tpgADC", &tpgADC_, "tpgADC/I");
0053 tree_->Branch("tpgGeV", &tpgGeV_, "tpgGeV/F");
0054 tree_->Branch("ttf", &ttf_, "ttf/I");
0055 tree_->Branch("fg", &fg_, "fg/I");
0056 for (unsigned int i = 0; i < 2; ++i) {
0057 ecal_et_[i] = fs->make<TH1I>(ecal_parts_[i].c_str(), "Et", 255, 0, 255);
0058 char title[30];
0059 sprintf(title, "%s_ttf", ecal_parts_[i].c_str());
0060 ecal_tt_[i] = fs->make<TH1I>(title, "TTF", 10, 0, 10);
0061 sprintf(title, "%s_fgvb", ecal_parts_[i].c_str());
0062 ecal_fgvb_[i] = fs->make<TH1I>(title, "FGVB", 10, 0, 10);
0063 }
0064
0065 if (recHits_) {
0066 hTPvsRechit_ = fs->make<TH2F>("TP_vs_RecHit", "TP vs rechit", 256, -1, 255, 255, 0, 255);
0067 hTPoverRechit_ = fs->make<TH1F>("TP_over_RecHit", "TP over rechit", 500, 0, 4);
0068 endcapGeomToken_ = esConsumes<CaloSubdetectorGeometry, EcalEndcapGeometryRecord>(edm::ESInputTag("", "EcalEndcap"));
0069 barrelGeomToken_ = esConsumes<CaloSubdetectorGeometry, EcalBarrelGeometryRecord>(edm::ESInputTag("", "EcalBarrel"));
0070 eTTmapToken_ = esConsumes<EcalTrigTowerConstituentsMap, IdealGeometryRecord>();
0071 }
0072 }
0073
0074
0075
0076
0077
0078
0079 void EcalTrigPrimAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0080
0081 const auto &tp = iEvent.get(tpToken_);
0082 for (unsigned int i = 0; i < tp.size(); i++) {
0083 const EcalTriggerPrimitiveDigi &d = tp[i];
0084 int subdet = d.id().subDet() - 1;
0085 if (subdet == 0) {
0086 ecal_et_[subdet]->Fill(d.compressedEt());
0087 } else {
0088 if (d.id().ietaAbs() == 27 || d.id().ietaAbs() == 28) {
0089 if (i % 2)
0090 ecal_et_[subdet]->Fill(d.compressedEt() * 2.);
0091 } else
0092 ecal_et_[subdet]->Fill(d.compressedEt());
0093 }
0094 ecal_tt_[subdet]->Fill(d.ttFlag());
0095 ecal_fgvb_[subdet]->Fill(d.fineGrain());
0096 }
0097 if (!recHits_)
0098 return;
0099
0100
0101 const EcalRecHitCollection &rechit_EB_col = iEvent.get(ebToken_);
0102 const EcalRecHitCollection &rechit_EE_col = iEvent.get(eeToken_);
0103
0104 const auto &theEndcapGeometry = iSetup.getData(endcapGeomToken_);
0105 const auto &theBarrelGeometry = iSetup.getData(barrelGeomToken_);
0106 const auto &eTTmap = iSetup.getData(eTTmapToken_);
0107
0108 std::map<EcalTrigTowerDetId, float> mapTow_Et;
0109
0110 for (unsigned int i = 0; i < rechit_EB_col.size(); i++) {
0111 const EBDetId &myid1 = rechit_EB_col[i].id();
0112 EcalTrigTowerDetId towid1 = myid1.tower();
0113 float theta = theBarrelGeometry.getGeometry(myid1)->getPosition().theta();
0114 float Etsum = rechit_EB_col[i].energy() * sin(theta);
0115 bool test_alreadyin = false;
0116 std::map<EcalTrigTowerDetId, float>::iterator ittest = mapTow_Et.find(towid1);
0117 if (ittest != mapTow_Et.end())
0118 test_alreadyin = true;
0119 if (test_alreadyin)
0120 continue;
0121 unsigned int j = i + 1;
0122 bool loopend = false;
0123 unsigned int count = 0;
0124 while (j < rechit_EB_col.size() && !loopend) {
0125 count++;
0126 const EBDetId &myid2 = rechit_EB_col[j].id();
0127 EcalTrigTowerDetId towid2 = myid2.tower();
0128 if (towid1 == towid2) {
0129 float theta = theBarrelGeometry.getGeometry(myid2)->getPosition().theta();
0130 Etsum += rechit_EB_col[j].energy() * sin(theta);
0131 }
0132 j++;
0133 if (count > 1800)
0134 loopend = true;
0135 }
0136 mapTow_Et.insert(std::pair<EcalTrigTowerDetId, float>(towid1, Etsum));
0137 }
0138
0139 for (unsigned int i = 0; i < rechit_EE_col.size(); i++) {
0140 const EEDetId &myid1 = rechit_EE_col[i].id();
0141 EcalTrigTowerDetId towid1 = eTTmap.towerOf(myid1);
0142 float theta = theEndcapGeometry.getGeometry(myid1)->getPosition().theta();
0143 float Etsum = rechit_EE_col[i].energy() * sin(theta);
0144 bool test_alreadyin = false;
0145 std::map<EcalTrigTowerDetId, float>::iterator ittest = mapTow_Et.find(towid1);
0146 if (ittest != mapTow_Et.end())
0147 test_alreadyin = true;
0148 if (test_alreadyin)
0149 continue;
0150 unsigned int j = i + 1;
0151 bool loopend = false;
0152 unsigned int count = 0;
0153 while (j < rechit_EE_col.size() && !loopend) {
0154 const EEDetId &myid2 = rechit_EE_col[j].id();
0155 EcalTrigTowerDetId towid2 = eTTmap.towerOf(myid2);
0156 if (towid1 == towid2) {
0157 float theta = theEndcapGeometry.getGeometry(myid2)->getPosition().theta();
0158 Etsum += rechit_EE_col[j].energy() * sin(theta);
0159 }
0160
0161 j++;
0162 if (count > 500)
0163 loopend = true;
0164 }
0165
0166 mapTow_Et.insert(std::pair<EcalTrigTowerDetId, float>(towid1, Etsum));
0167 }
0168
0169 EcalTPGScale ecalScale(tokens_, iSetup);
0170 for (unsigned int i = 0; i < tp.size(); i++) {
0171 const EcalTriggerPrimitiveDigi &d = tp[i];
0172 const EcalTrigTowerDetId TPtowid = d.id();
0173 std::map<EcalTrigTowerDetId, float>::iterator it = mapTow_Et.find(TPtowid);
0174 float Et = ecalScale.getTPGInGeV(d.compressedEt(), TPtowid);
0175 if (d.id().ietaAbs() == 27 || d.id().ietaAbs() == 28)
0176 Et *= 2;
0177 iphi_ = TPtowid.iphi();
0178 ieta_ = TPtowid.ieta();
0179 tpgADC_ = d.compressedEt();
0180 tpgGeV_ = Et;
0181 ttf_ = d.ttFlag();
0182 fg_ = d.fineGrain();
0183 if (it != mapTow_Et.end()) {
0184 hTPvsRechit_->Fill(it->second, Et);
0185 hTPoverRechit_->Fill(Et / it->second);
0186 eRec_ = it->second;
0187 }
0188 tree_->Fill();
0189 }
0190 }
