File indexing completed on 2024-09-07 04:34:56
0001
0002
0003
0004
0005
0006
0007 #include "CalibMuon/DTCalibration/plugins/DTNoiseComputation.h"
0008 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0009
0010
0011 #include "FWCore/Framework/interface/IOVSyncValue.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016
0017
0018 #include "Geometry/DTGeometry/interface/DTLayer.h"
0019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0020 #include "Geometry/DTGeometry/interface/DTTopology.h"
0021 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0022
0023
0024 #include "DataFormats/DTDigi/interface/DTDigi.h"
0025 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0026 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
0027 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0028
0029 #include "TH1F.h"
0030 #include "TH2F.h"
0031 #include "TFile.h"
0032 #include "TF1.h"
0033 #include "TProfile.h"
0034 #include "TPostScript.h"
0035 #include "TCanvas.h"
0036 #include "TLegend.h"
0037
0038 using namespace edm;
0039 using namespace std;
0040
0041 DTNoiseComputation::DTNoiseComputation(const edm::ParameterSet &ps) : dtGeomToken_(esConsumes()) {
0042 cout << "[DTNoiseComputation]: Constructor" << endl;
0043
0044
0045 debug = ps.getUntrackedParameter<bool>("debug");
0046
0047
0048 fastAnalysis = ps.getUntrackedParameter<bool>("fastAnalysis", true);
0049
0050
0051 string rootFileName = ps.getUntrackedParameter<string>("rootFileName");
0052 theFile = new TFile(rootFileName.c_str(), "READ");
0053
0054
0055 string newRootFileName = ps.getUntrackedParameter<string>("newRootFileName");
0056 theNewFile = new TFile(newRootFileName.c_str(), "RECREATE");
0057
0058
0059 MaxEvents = ps.getUntrackedParameter<int>("MaxEvents");
0060 }
0061
0062 void DTNoiseComputation::beginRun(const edm::Run &, const EventSetup &setup) {
0063
0064 dtGeom = setup.getHandle(dtGeomToken_);
0065
0066 static int count = 0;
0067
0068 if (count == 0) {
0069 string CheckHistoName;
0070
0071 TH1F *hOccHisto;
0072 TH1F *hAverageNoiseHisto;
0073 TH1F *hAverageNoiseIntegratedHisto;
0074 TH1F *hAverageNoiseHistoPerCh;
0075 TH1F *hAverageNoiseIntegratedHistoPerCh;
0076 TH2F *hEvtHisto;
0077 string HistoName;
0078 string Histo2Name;
0079 string AverageNoiseName;
0080 string AverageNoiseIntegratedName;
0081 string AverageNoiseNamePerCh;
0082 string AverageNoiseIntegratedNamePerCh;
0083 TH1F *hnoisyC;
0084 TH1F *hsomeHowNoisyC;
0085
0086
0087 vector<const DTChamber *>::const_iterator ch_it = dtGeom->chambers().begin();
0088 vector<const DTChamber *>::const_iterator ch_end = dtGeom->chambers().end();
0089
0090 for (; ch_it != ch_end; ++ch_it) {
0091 DTChamberId ch = (*ch_it)->id();
0092 vector<const DTSuperLayer *>::const_iterator sl_it = (*ch_it)->superLayers().begin();
0093 vector<const DTSuperLayer *>::const_iterator sl_end = (*ch_it)->superLayers().end();
0094
0095 for (; sl_it != sl_end; ++sl_it) {
0096
0097 vector<const DTLayer *>::const_iterator l_it = (*sl_it)->layers().begin();
0098 vector<const DTLayer *>::const_iterator l_end = (*sl_it)->layers().end();
0099
0100 for (; l_it != l_end; ++l_it) {
0101 DTLayerId dtLId = (*l_it)->id();
0102
0103
0104 theFile->cd();
0105 CheckHistoName = "DigiOccupancy_" + getLayerName(dtLId);
0106 TH1F *hCheckHisto = (TH1F *)theFile->Get(CheckHistoName.c_str());
0107 if (hCheckHisto) {
0108 delete hCheckHisto;
0109 string wheel = std::to_string(ch.wheel());
0110 string station = std::to_string(ch.station());
0111
0112 if (someHowNoisyC.find(make_pair(ch.wheel(), ch.station())) == someHowNoisyC.end()) {
0113 TString histoName_someHowNoisy = "somehowNoisyCell_W" + wheel + "_St" + station;
0114 hsomeHowNoisyC =
0115 new TH1F(histoName_someHowNoisy, histoName_someHowNoisy, getMaxNumBins(ch), 1, getMaxNumBins(ch) + 1);
0116 someHowNoisyC[make_pair(ch.wheel(), ch.station())] = hsomeHowNoisyC;
0117 }
0118
0119 if (noisyC.find(make_pair(ch.wheel(), ch.station())) == noisyC.end()) {
0120 TString histoName_noisy = "noisyCell_W" + wheel + "_St" + station;
0121 hnoisyC = new TH1F(histoName_noisy, histoName_noisy, getMaxNumBins(ch), 1, getMaxNumBins(ch) + 1);
0122 noisyC[make_pair(ch.wheel(), ch.station())] = hnoisyC;
0123 }
0124
0125
0126 if (AvNoisePerSuperLayer.find(dtLId.superlayerId()) == AvNoisePerSuperLayer.end()) {
0127 AverageNoiseName = "AverageNoise_" + getSuperLayerName(dtLId.superlayerId());
0128 hAverageNoiseHisto = new TH1F(AverageNoiseName.c_str(), AverageNoiseName.c_str(), 200, 0, 10000);
0129 AverageNoiseIntegratedName = "AverageNoiseIntegrated_" + getSuperLayerName(dtLId.superlayerId());
0130 hAverageNoiseIntegratedHisto =
0131 new TH1F(AverageNoiseIntegratedName.c_str(), AverageNoiseIntegratedName.c_str(), 200, 0, 10000);
0132 AvNoisePerSuperLayer[dtLId.superlayerId()] = hAverageNoiseHisto;
0133 AvNoiseIntegratedPerSuperLayer[dtLId.superlayerId()] = hAverageNoiseIntegratedHisto;
0134 if (debug) {
0135 cout << " New Average Noise Histo per SuperLayer : " << hAverageNoiseHisto->GetName() << endl;
0136 cout << " New Average Noise Integrated Histo per SuperLayer : " << hAverageNoiseHisto->GetName()
0137 << endl;
0138 }
0139 }
0140 if (AvNoisePerChamber.find(dtLId.superlayerId().chamberId()) == AvNoisePerChamber.end()) {
0141 AverageNoiseNamePerCh = "AverageNoise_" + getChamberName(dtLId);
0142 hAverageNoiseHistoPerCh =
0143 new TH1F(AverageNoiseNamePerCh.c_str(), AverageNoiseNamePerCh.c_str(), 200, 0, 10000);
0144 AverageNoiseIntegratedNamePerCh = "AverageNoiseIntegrated_" + getChamberName(dtLId);
0145 hAverageNoiseIntegratedHistoPerCh = new TH1F(
0146 AverageNoiseIntegratedNamePerCh.c_str(), AverageNoiseIntegratedNamePerCh.c_str(), 200, 0, 10000);
0147 AvNoisePerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseHistoPerCh;
0148 AvNoiseIntegratedPerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseIntegratedHistoPerCh;
0149 if (debug)
0150 cout << " New Average Noise Histo per chamber : " << hAverageNoiseHistoPerCh->GetName() << endl;
0151 }
0152
0153 HistoName = "DigiOccupancy_" + getLayerName(dtLId);
0154 theFile->cd();
0155 hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0156 int numBin = hOccHisto->GetXaxis()->GetNbins();
0157 for (int bin = 1; bin <= numBin; bin++) {
0158 DTWireId wireID(dtLId, bin);
0159 theAverageNoise[wireID] = hOccHisto->GetBinContent(bin);
0160 if (theAverageNoise[wireID] != 0) {
0161 AvNoisePerSuperLayer[dtLId.superlayerId()]->Fill(theAverageNoise[wireID]);
0162 AvNoisePerChamber[dtLId.superlayerId().chamberId()]->Fill(theAverageNoise[wireID]);
0163 }
0164 }
0165
0166
0167 double numCell = 0;
0168 double AvNoise = 0;
0169 HistoName = "DigiOccupancy_" + getLayerName(dtLId);
0170 theFile->cd();
0171 hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0172 numBin = hOccHisto->GetXaxis()->GetNbins();
0173 for (int bin = 1; bin <= numBin; bin++) {
0174 DTWireId wireID(dtLId, bin);
0175 theAverageNoise[wireID] = hOccHisto->GetBinContent(bin);
0176 if (hOccHisto->GetBinContent(bin) < 100) {
0177 numCell++;
0178 AvNoise += hOccHisto->GetBinContent(bin);
0179 }
0180 if (hOccHisto->GetBinContent(bin) > 100 && hOccHisto->GetBinContent(bin) < 500) {
0181 someHowNoisyC[make_pair(ch.wheel(), ch.station())]->Fill(bin);
0182 cout << "filling somehow noisy cell" << endl;
0183 }
0184 if (hOccHisto->GetBinContent(bin) > 500) {
0185 noisyC[make_pair(ch.wheel(), ch.station())]->Fill(bin);
0186 cout << "filling noisy cell" << endl;
0187 }
0188 }
0189 AvNoise = AvNoise / numCell;
0190 cout << "theAverageNoise for layer " << getLayerName(dtLId) << " is : " << AvNoise << endl;
0191
0192
0193 int updates = MaxEvents / 1000;
0194 for (int evt = 0; evt < updates; evt++) {
0195 Histo2Name = "DigiPerWirePerEvent_" + getLayerName(dtLId) + "_" + std::to_string(evt);
0196 theFile->cd();
0197 hEvtHisto = (TH2F *)theFile->Get(Histo2Name.c_str());
0198 if (hEvtHisto) {
0199 if (debug)
0200 cout << " New Histo with the number of events per evt per wire: " << hEvtHisto->GetName() << endl;
0201 theEvtMap[dtLId].push_back(hEvtHisto);
0202 }
0203 }
0204
0205 }
0206 }
0207 }
0208 }
0209
0210 count++;
0211 }
0212 }
0213
0214 void DTNoiseComputation::endJob() {
0215 cout << "[DTNoiseComputation] endjob called!" << endl;
0216 TH1F *hEvtDistance = nullptr;
0217 TF1 *ExpoFit = new TF1("ExpoFit", "expo", 0.5, 1000.5);
0218 ExpoFit->SetMarkerColor();
0219 TProfile *theNoiseHisto = new TProfile("theNoiseHisto", "Time Constant versus Average Noise", 100000, 0, 100000);
0220
0221
0222 for (map<DTLayerId, vector<TH2F *> >::const_iterator lHisto = theEvtMap.begin(); lHisto != theEvtMap.end();
0223 ++lHisto) {
0224 for (int bin = 1; bin < (*lHisto).second[0]->GetYaxis()->GetNbins(); bin++) {
0225 int distanceEvt = 1;
0226 DTWireId wire((*lHisto).first, bin);
0227 for (int i = 0; i < int((*lHisto).second.size()); i++) {
0228 for (int evt = 1; evt <= (*lHisto).second[i]->GetXaxis()->GetNbins(); evt++) {
0229 if ((*lHisto).second[i]->GetBinContent(evt, bin) == 0)
0230 distanceEvt++;
0231 else {
0232 if (toDel.find(wire) == toDel.end()) {
0233 toDel[wire] = false;
0234 string Histo = "EvtDistancePerWire_" + getLayerName((*lHisto).first) + "_" + std::to_string(bin);
0235 hEvtDistance = new TH1F(Histo.c_str(), Histo.c_str(), 50000, 0.5, 50000.5);
0236 }
0237 hEvtDistance->Fill(distanceEvt);
0238 distanceEvt = 1;
0239 }
0240 }
0241 }
0242 if (toDel.find(wire) != toDel.end()) {
0243 theHistoEvtDistancePerWire[wire] = hEvtDistance;
0244 theNewFile->cd();
0245 theHistoEvtDistancePerWire[wire]->Fit("ExpoFit", "R");
0246 TF1 *funct = theHistoEvtDistancePerWire[wire]->GetFunction("ExpoFit");
0247 double par0 = funct->GetParameter(0);
0248 double par1 = funct->GetParameter(1);
0249 cout << "par0: " << par0 << " par1: " << par1 << endl;
0250 double chi2rid = funct->GetChisquare() / funct->GetNDF();
0251 if (chi2rid > 10)
0252 theTimeConstant[wire] = 1;
0253 else
0254 theTimeConstant[wire] = par1;
0255 toDel[wire] = true;
0256 theHistoEvtDistancePerWire[wire]->Write();
0257 delete hEvtDistance;
0258 }
0259 }
0260 }
0261
0262 if (!fastAnalysis) {
0263
0264 for (map<DTWireId, double>::const_iterator AvNoise = theAverageNoise.begin(); AvNoise != theAverageNoise.end();
0265 ++AvNoise) {
0266 DTWireId wire = (*AvNoise).first;
0267 theNoiseHisto->Fill((*AvNoise).second, theTimeConstant[wire]);
0268 cout << "Layer: " << getLayerName(wire.layerId()) << " wire: " << wire.wire() << endl;
0269 cout << "The Average noise: " << (*AvNoise).second << endl;
0270 cout << "The time constant: " << theTimeConstant[wire] << endl;
0271 }
0272 theNewFile->cd();
0273 theNoiseHisto->Write();
0274 }
0275
0276
0277 int numBin;
0278 double integratedNoise, bin, halfBin, maxBin;
0279 for (map<DTSuperLayerId, TH1F *>::const_iterator AvNoiseHisto = AvNoisePerSuperLayer.begin();
0280 AvNoiseHisto != AvNoisePerSuperLayer.end();
0281 ++AvNoiseHisto) {
0282 integratedNoise = 0;
0283 numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
0284 maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
0285 bin = double(maxBin / numBin);
0286 halfBin = double(bin / 2);
0287 theNewFile->cd();
0288 (*AvNoiseHisto).second->Write();
0289 for (int i = 1; i < numBin; i++) {
0290 integratedNoise += (*AvNoiseHisto).second->GetBinContent(i);
0291 AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Fill(halfBin, integratedNoise);
0292 halfBin += bin;
0293 }
0294 theNewFile->cd();
0295 AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Write();
0296 }
0297
0298 for (map<DTChamberId, TH1F *>::const_iterator AvNoiseHisto = AvNoisePerChamber.begin();
0299 AvNoiseHisto != AvNoisePerChamber.end();
0300 ++AvNoiseHisto) {
0301 integratedNoise = 0;
0302 numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
0303 maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
0304 bin = maxBin / numBin;
0305 halfBin = bin / 2;
0306 theNewFile->cd();
0307 (*AvNoiseHisto).second->Write();
0308 for (int i = 1; i < numBin; i++) {
0309 integratedNoise += (*AvNoiseHisto).second->GetBinContent(i);
0310 AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Fill(halfBin, integratedNoise);
0311 halfBin += bin;
0312 }
0313 theNewFile->cd();
0314 AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Write();
0315 }
0316
0317
0318 bool histo = false;
0319 vector<const DTChamber *>::const_iterator chamber_it = dtGeom->chambers().begin();
0320 vector<const DTChamber *>::const_iterator chamber_end = dtGeom->chambers().end();
0321
0322 for (; chamber_it != chamber_end; ++chamber_it) {
0323 vector<const DTSuperLayer *>::const_iterator sl_it = (*chamber_it)->superLayers().begin();
0324 vector<const DTSuperLayer *>::const_iterator sl_end = (*chamber_it)->superLayers().end();
0325
0326 for (; sl_it != sl_end; ++sl_it) {
0327 DTSuperLayerId sl = (*sl_it)->id();
0328 vector<const DTLayer *>::const_iterator l_it = (*sl_it)->layers().begin();
0329 vector<const DTLayer *>::const_iterator l_end = (*sl_it)->layers().end();
0330
0331 string canvasName = "c" + getSuperLayerName(sl);
0332 TCanvas c1(canvasName.c_str(), canvasName.c_str(), 600, 780);
0333 TLegend *leg = new TLegend(0.5, 0.6, 0.7, 0.8);
0334 for (; l_it != l_end; ++l_it) {
0335 DTLayerId layerId = (*l_it)->id();
0336 string HistoName = "DigiOccupancy_" + getLayerName(layerId);
0337 theFile->cd();
0338 TH1F *hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0339 if (hOccHisto) {
0340 string TitleHisto = "AverageNoise_" + getSuperLayerName(sl);
0341 cout << "TitleHisto : " << TitleHisto << endl;
0342 hOccHisto->SetTitle(TitleHisto.c_str());
0343 string legendHisto = "layer " + std::to_string(layerId.layer());
0344 leg->AddEntry(hOccHisto, legendHisto.c_str(), "L");
0345 hOccHisto->SetMaximum(getYMaximum(sl));
0346 histo = true;
0347 if (layerId.layer() == 1)
0348 hOccHisto->Draw();
0349 else
0350 hOccHisto->Draw("same");
0351 hOccHisto->SetLineColor(layerId.layer());
0352 }
0353 }
0354 if (histo) {
0355 leg->Draw("same");
0356 theNewFile->cd();
0357 c1.Write();
0358 }
0359 }
0360 histo = false;
0361 }
0362
0363
0364 for (map<pair<int, int>, TH1F *>::const_iterator nCell = noisyC.begin(); nCell != noisyC.end(); ++nCell) {
0365 theNewFile->cd();
0366 (*nCell).second->Write();
0367 }
0368 for (map<pair<int, int>, TH1F *>::const_iterator somehownCell = someHowNoisyC.begin();
0369 somehownCell != someHowNoisyC.end();
0370 ++somehownCell) {
0371 theNewFile->cd();
0372 (*somehownCell).second->Write();
0373 }
0374 }
0375
0376 DTNoiseComputation::~DTNoiseComputation() {
0377 theFile->Close();
0378 theNewFile->Close();
0379 }
0380
0381 string DTNoiseComputation::getLayerName(const DTLayerId &lId) const {
0382 const DTSuperLayerId dtSLId = lId.superlayerId();
0383 const DTChamberId dtChId = dtSLId.chamberId();
0384 string layerName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0385 std::to_string(dtChId.sector()) + "_SL" + std::to_string(dtSLId.superlayer()) + "_L" +
0386 std::to_string(lId.layer());
0387
0388 return layerName;
0389 }
0390
0391 string DTNoiseComputation::getSuperLayerName(const DTSuperLayerId &dtSLId) const {
0392 const DTChamberId dtChId = dtSLId.chamberId();
0393
0394 string superLayerName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0395 std::to_string(dtChId.sector()) + "_SL" + std::to_string(dtSLId.superlayer());
0396
0397 return superLayerName;
0398 }
0399
0400 string DTNoiseComputation::getChamberName(const DTLayerId &lId) const {
0401 const DTSuperLayerId dtSLId = lId.superlayerId();
0402 const DTChamberId dtChId = dtSLId.chamberId();
0403 string chamberName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0404 std::to_string(dtChId.sector());
0405
0406 return chamberName;
0407 }
0408
0409 int DTNoiseComputation::getMaxNumBins(const DTChamberId &chId) const {
0410 int maximum = 0;
0411
0412 for (int SL = 1; SL <= 3; SL++) {
0413 if (!(chId.station() == 4 && SL == 2)) {
0414 for (int L = 1; L <= 4; L++) {
0415 DTLayerId layerId = DTLayerId(chId, SL, L);
0416 string HistoName = "DigiOccupancy_" + getLayerName(layerId);
0417 theFile->cd();
0418 TH1F *hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0419 if (hOccHisto) {
0420 if (hOccHisto->GetXaxis()->GetXmax() > maximum)
0421 maximum = hOccHisto->GetXaxis()->GetNbins();
0422 }
0423 }
0424 }
0425 }
0426 return maximum;
0427 }
0428
0429 double DTNoiseComputation::getYMaximum(const DTSuperLayerId &slId) const {
0430 double maximum = 0;
0431 double dummy = pow(10., 10.);
0432
0433 for (int L = 1; L <= 4; L++) {
0434 DTLayerId layerId = DTLayerId(slId, L);
0435 string HistoName = "DigiOccupancy_" + getLayerName(layerId);
0436 theFile->cd();
0437 TH1F *hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0438 if (hOccHisto) {
0439 if (hOccHisto->GetMaximum(dummy) > maximum)
0440 maximum = hOccHisto->GetMaximum(dummy);
0441 }
0442 }
0443 return maximum;
0444 }