File indexing completed on 2023-03-17 10:44:49
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include "CaloOnlineTools/EcalTools/plugins/EcalMipGraphs.h"
0020 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0021 #include "TCanvas.h"
0022 #include <utility>
0023
0024 using namespace edm;
0025 using namespace std;
0026
0027
0028
0029
0030
0031
0032
0033
0034 float EcalMipGraphs::gainRatio[3] = {1., 2., 12.};
0035 edm::Service<TFileService> EcalMipGraphs::fileService;
0036
0037
0038
0039
0040 EcalMipGraphs::EcalMipGraphs(const edm::ParameterSet& iConfig)
0041 : EBRecHitCollection_(iConfig.getParameter<edm::InputTag>("EcalRecHitCollectionEB")),
0042 EERecHitCollection_(iConfig.getParameter<edm::InputTag>("EcalRecHitCollectionEE")),
0043 EBDigis_(iConfig.getParameter<edm::InputTag>("EBDigiCollection")),
0044 EEDigis_(iConfig.getParameter<edm::InputTag>("EEDigiCollection")),
0045 headerProducer_(iConfig.getParameter<edm::InputTag>("headerProducer")),
0046 rawDataToken_(consumes<EcalRawDataCollection>(headerProducer_)),
0047 ebRecHitToken_(consumes<EcalRecHitCollection>(EBRecHitCollection_)),
0048 eeRecHitToken_(consumes<EcalRecHitCollection>(EERecHitCollection_)),
0049 ebDigiToken_(consumes<EBDigiCollection>(EBDigis_)),
0050 eeDigiToken_(consumes<EEDigiCollection>(EEDigis_)),
0051 ecalMappingToken_(esConsumes<edm::Transition::BeginRun>()),
0052 topologyToken_(esConsumes()),
0053 runNum_(-1),
0054 side_(iConfig.getUntrackedParameter<int>("side", 3)),
0055 threshold_(iConfig.getUntrackedParameter<double>("amplitudeThreshold", 12.0)),
0056 minTimingAmp_(iConfig.getUntrackedParameter<double>("minimumTimingAmplitude", 0.100)) {
0057 vector<int> listDefaults;
0058 listDefaults.push_back(-1);
0059
0060 maskedChannels_ = iConfig.getUntrackedParameter<vector<int> >("maskedChannels", listDefaults);
0061 maskedFEDs_ = iConfig.getUntrackedParameter<vector<int> >("maskedFEDs", listDefaults);
0062 seedCrys_ = iConfig.getUntrackedParameter<vector<int> >("seedCrys", vector<int>());
0063
0064 vector<string> defaultMaskedEBs;
0065 defaultMaskedEBs.push_back("none");
0066 maskedEBs_ = iConfig.getUntrackedParameter<vector<string> >("maskedEBs", defaultMaskedEBs);
0067
0068 fedMap_ = new EcalFedMap();
0069
0070 string title1 = "Jitter for all FEDs";
0071 string name1 = "JitterAllFEDs";
0072 allFedsTimingHist_ = fileService->make<TH1F>(name1.c_str(), title1.c_str(), 150, -7, 7);
0073
0074
0075 if (maskedFEDs_[0] == -1) {
0076
0077 if (maskedEBs_[0] != "none") {
0078 maskedFEDs_.clear();
0079 for (vector<string>::const_iterator ebItr = maskedEBs_.begin(); ebItr != maskedEBs_.end(); ++ebItr) {
0080 maskedFEDs_.push_back(fedMap_->getFedFromSlice(*ebItr));
0081 }
0082 }
0083 }
0084
0085 for (int i = 0; i < 10; i++)
0086 abscissa[i] = i;
0087
0088 naiveEvtNum_ = 0;
0089 }
0090
0091 EcalMipGraphs::~EcalMipGraphs() {}
0092
0093
0094
0095
0096
0097
0098 void EcalMipGraphs::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0099
0100
0101 edm::Handle<EcalRawDataCollection> DCCHeaders;
0102 iEvent.getByToken(rawDataToken_, DCCHeaders);
0103
0104 for (EcalRawDataCollection::const_iterator headerItr = DCCHeaders->begin(); headerItr != DCCHeaders->end();
0105 ++headerItr) {
0106 FEDsAndDCCHeaders_[headerItr->id() + 600] = *headerItr;
0107 }
0108
0109 int ievt = iEvent.id().event();
0110 naiveEvtNum_++;
0111
0112 if (runNum_ == -1) {
0113 runNum_ = iEvent.id().run();
0114 canvasNames_ = fileService->make<TTree>("canvasNames", "Names of written canvases");
0115 names = new std::vector<string>();
0116 canvasNames_->Branch("names", "vector<string>", &names);
0117 }
0118
0119
0120 listEBChannels.clear();
0121 listEEChannels.clear();
0122 caloTopo_ = &iSetup.getData(topologyToken_);
0123
0124 Handle<EcalRecHitCollection> EBhits;
0125 Handle<EcalRecHitCollection> EEhits;
0126 iEvent.getByToken(ebRecHitToken_, EBhits);
0127 iEvent.getByToken(eeRecHitToken_, EEhits);
0128
0129 iEvent.getByToken(ebDigiToken_, EBdigisHandle);
0130 iEvent.getByToken(eeDigiToken_, EEdigisHandle);
0131
0132 selectHits(EBhits, ievt);
0133 selectHits(EEhits, ievt);
0134 }
0135
0136 TGraph* EcalMipGraphs::selectDigi(DetId thisDet, int ievt) {
0137 int emptyY[10];
0138 for (int i = 0; i < 10; i++)
0139 emptyY[i] = 0;
0140 TGraph* emptyGraph = fileService->make<TGraph>(10, abscissa, emptyY);
0141 emptyGraph->SetTitle("NOT ECAL");
0142
0143
0144 if (thisDet.det() != DetId::Ecal)
0145 return emptyGraph;
0146
0147 emptyGraph->SetTitle("NO DIGIS");
0148
0149 EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(thisDet);
0150 int FEDid = 600 + elecId.dccId();
0151 bool isBarrel = true;
0152 if (FEDid < 610 || FEDid > 645)
0153 isBarrel = false;
0154 int cryIndex = isBarrel ? ((EBDetId)thisDet).hashedIndex() : getEEIndex(elecId);
0155 int ic = isBarrel ? ((EBDetId)thisDet).ic() : cryIndex;
0156
0157 string sliceName = fedMap_->getSliceFromFed(FEDid);
0158 EcalDataFrame df;
0159 if (isBarrel) {
0160 EBDigiCollection::const_iterator digiItr = EBdigisHandle->begin();
0161 while (digiItr != EBdigisHandle->end() && ((*digiItr).id() != (EBDetId)thisDet)) {
0162 ++digiItr;
0163 }
0164 if (digiItr == EBdigisHandle->end()) {
0165
0166
0167 return emptyGraph;
0168 } else
0169 df = *digiItr;
0170 } else {
0171 EEDigiCollection::const_iterator digiItr = EEdigisHandle->begin();
0172 while (digiItr != EEdigisHandle->end() && ((*digiItr).id() != (EEDetId)thisDet)) {
0173 ++digiItr;
0174 }
0175 if (digiItr == EEdigisHandle->end()) {
0176
0177
0178 return emptyGraph;
0179 } else
0180 df = *digiItr;
0181 }
0182
0183 int gainId = FEDsAndDCCHeaders_[FEDid].getMgpaGain();
0184 int gainHuman;
0185 if (gainId == 1)
0186 gainHuman = 12;
0187 else if (gainId == 2)
0188 gainHuman = 6;
0189 else if (gainId == 3)
0190 gainHuman = 1;
0191 else
0192 gainHuman = -1;
0193
0194 double pedestal = 200;
0195
0196 emptyGraph->SetTitle("FIRST TWO SAMPLES NOT GAIN12");
0197 if (df.sample(0).gainId() != 1 || df.sample(1).gainId() != 1)
0198 return emptyGraph;
0199 else {
0200 ordinate[0] = df.sample(0).adc();
0201 ordinate[1] = df.sample(1).adc();
0202 pedestal = (double)(ordinate[0] + ordinate[1]) / (double)2;
0203 }
0204
0205 for (int i = 0; i < df.size(); ++i) {
0206 if (df.sample(i).gainId() != 0)
0207 ordinate[i] = (int)(pedestal + (df.sample(i).adc() - pedestal) * gainRatio[df.sample(i).gainId() - 1]);
0208 else
0209 ordinate[i] = 49152;
0210 }
0211
0212 TGraph* oneGraph = fileService->make<TGraph>(10, abscissa, ordinate);
0213 string name = "Graph_ev" + intToString(naiveEvtNum_) + "_ic" + intToString(ic) + "_FED" + intToString(FEDid);
0214 string gainString = (gainId == 1) ? "Free" : intToString(gainHuman);
0215 string title = "Event" + intToString(naiveEvtNum_) + "_lv1a" + intToString(ievt) + "_ic" + intToString(ic) + "_" +
0216 sliceName + "_gain" + gainString;
0217
0218 float energy = 0;
0219 map<int, float>::const_iterator itr;
0220 itr = crysAndAmplitudesMap_.find(cryIndex);
0221 if (itr != crysAndAmplitudesMap_.end()) {
0222
0223 energy = (float)itr->second;
0224 }
0225
0226
0227
0228 title += "_Energy" + floatToString(round(energy * 1000));
0229
0230 oneGraph->SetTitle(title.c_str());
0231 oneGraph->SetName(name.c_str());
0232 oneGraph->GetXaxis()->SetTitle("sample");
0233 oneGraph->GetYaxis()->SetTitle("ADC");
0234 return oneGraph;
0235 }
0236
0237 void EcalMipGraphs::selectHits(Handle<EcalRecHitCollection> hits, int ievt) {
0238 for (EcalRecHitCollection::const_iterator hitItr = hits->begin(); hitItr != hits->end(); ++hitItr) {
0239 EcalRecHit hit = (*hitItr);
0240 DetId det = hit.id();
0241 EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(det);
0242 int FEDid = 600 + elecId.dccId();
0243 bool isBarrel = true;
0244 if (FEDid < 610 || FEDid > 645)
0245 isBarrel = false;
0246 int cryIndex = isBarrel ? ((EBDetId)det).hashedIndex() : ((EEDetId)det).hashedIndex();
0247 int ic = isBarrel ? ((EBDetId)det).ic() : getEEIndex(elecId);
0248
0249 float ampli = hit.energy();
0250
0251 vector<int>::iterator result;
0252 result = find(maskedFEDs_.begin(), maskedFEDs_.end(), FEDid);
0253 if (result != maskedFEDs_.end()) {
0254
0255 continue;
0256 }
0257 result = find(maskedChannels_.begin(), maskedChannels_.end(), cryIndex);
0258 if (result != maskedChannels_.end()) {
0259
0260 continue;
0261 }
0262 bool cryIsInList = false;
0263 result = find(seedCrys_.begin(), seedCrys_.end(), cryIndex);
0264 if (result != seedCrys_.end())
0265 cryIsInList = true;
0266
0267
0268
0269 if (cryIsInList || (seedCrys_.empty() && ampli > threshold_)) {
0270
0271 crysAndAmplitudesMap_[cryIndex] = ampli;
0272 string name = "Event" + intToString(naiveEvtNum_) + "_ic" + intToString(ic) + "_FED" + intToString(FEDid);
0273 string title = "Digis";
0274 string seed = "ic" + intToString(ic) + "_FED" + intToString(FEDid);
0275 int freq = 1;
0276 pair<map<string, int>::iterator, bool> pair = seedFrequencyMap_.insert(make_pair(seed, freq));
0277 if (!pair.second) {
0278 ++(pair.first->second);
0279 }
0280
0281 TCanvas can(name.c_str(), title.c_str(), 200, 50, 900, 900);
0282 can.Divide(side_, side_);
0283 TGraph* myGraph;
0284 int canvasNum = 1;
0285
0286 CaloNavigator<DetId> cursor = CaloNavigator<DetId>(det, caloTopo_->getSubdetectorTopology(det));
0287
0288 for (int j = side_ / 2; j >= -side_ / 2; --j) {
0289 for (int i = -side_ / 2; i <= side_ / 2; ++i) {
0290 cursor.home();
0291 cursor.offsetBy(i, j);
0292 can.cd(canvasNum);
0293 myGraph = selectDigi(*cursor, ievt);
0294 myGraph->Draw("A*");
0295 canvasNum++;
0296 }
0297 }
0298 can.Write();
0299 names->push_back(name);
0300 }
0301
0302 TH1F* timingHist = FEDsAndTimingHists_[FEDid];
0303 if (timingHist == nullptr) {
0304 initHists(FEDid);
0305 timingHist = FEDsAndTimingHists_[FEDid];
0306 }
0307 if (ampli > minTimingAmp_) {
0308 timingHist->Fill(hit.time());
0309 allFedsTimingHist_->Fill(hit.time());
0310 }
0311 }
0312 }
0313
0314 int EcalMipGraphs::getEEIndex(EcalElectronicsId elecId) {
0315 int FEDid = 600 + elecId.dccId();
0316 return 10000 * FEDid + 100 * elecId.towerId() + 5 * (elecId.stripId() - 1) + elecId.xtalId();
0317 }
0318
0319
0320 void EcalMipGraphs::initHists(int FED) {
0321 using namespace std;
0322
0323 string title1 = "Jitter for ";
0324 title1.append(fedMap_->getSliceFromFed(FED));
0325 string name1 = "JitterFED";
0326 name1.append(intToString(FED));
0327 TH1F* timingHist = fileService->make<TH1F>(name1.c_str(), title1.c_str(), 150, -7, 7);
0328 FEDsAndTimingHists_[FED] = timingHist;
0329 }
0330
0331
0332 void EcalMipGraphs::beginRun(edm::Run const&, edm::EventSetup const& c) {
0333 ecalElectronicsMap_ = &c.getData(ecalMappingToken_);
0334 }
0335
0336 void EcalMipGraphs::endRun(edm::Run const&, edm::EventSetup const& c) {}
0337
0338
0339 void EcalMipGraphs::endJob() {
0340 canvasNames_->Fill();
0341
0342 string frequencies = "";
0343 for (std::map<std::string, int>::const_iterator itr = seedFrequencyMap_.begin(); itr != seedFrequencyMap_.end();
0344 ++itr) {
0345 if (itr->second > 1) {
0346 frequencies += itr->first;
0347 frequencies += " Frequency: ";
0348 frequencies += intToString(itr->second);
0349 frequencies += "\n";
0350 }
0351 }
0352 LogWarning("EcalMipGraphs") << "Found seeds with frequency > 1: "
0353 << "\n\n"
0354 << frequencies;
0355
0356 std::string channels;
0357 for (std::vector<int>::const_iterator itr = maskedChannels_.begin(); itr != maskedChannels_.end(); ++itr) {
0358 channels += intToString(*itr);
0359 channels += ",";
0360 }
0361
0362 std::string feds;
0363 for (std::vector<int>::const_iterator itr = maskedFEDs_.begin(); itr != maskedFEDs_.end(); ++itr) {
0364 feds += intToString(*itr);
0365 feds += ",";
0366 }
0367
0368 LogWarning("EcalMipGraphs") << "Masked channels are: " << channels;
0369 LogWarning("EcalMipGraphs") << "Masked FEDs are: " << feds << " and that is all!";
0370 }
0371
0372 std::string EcalMipGraphs::intToString(int num) {
0373 using namespace std;
0374 ostringstream myStream;
0375 myStream << num << flush;
0376 return (myStream.str());
0377 }
0378
0379 std::string EcalMipGraphs::floatToString(float num) {
0380 using namespace std;
0381 ostringstream myStream;
0382 myStream << num << flush;
0383 return (myStream.str());
0384 }