File indexing completed on 2023-10-25 09:34:05
0001
0002
0003
0004
0005
0006 #include "CalibMuon/DTCalibration/plugins/DTT0CalibrationRMS.h"
0007 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0008
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0014 #include "Geometry/Records/interface/MuonNumberingRecord.h"
0015
0016 #include "CondFormats/DTObjects/interface/DTT0.h"
0017
0018 #include "TH1I.h"
0019 #include "TFile.h"
0020 #include "TKey.h"
0021
0022 using namespace std;
0023 using namespace edm;
0024
0025
0026
0027 DTT0CalibrationRMS::DTT0CalibrationRMS(const edm::ParameterSet& pset) : dtGeomToken_(esConsumes()) {
0028
0029 debug = pset.getUntrackedParameter<bool>("debug");
0030 if (debug)
0031 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS]Constructor called!";
0032
0033
0034 digiToken = consumes<DTDigiCollection>(pset.getUntrackedParameter<string>("digiLabel"));
0035
0036
0037 string rootFileName = pset.getUntrackedParameter<string>("rootFileName", "DTT0PerLayer.root");
0038 theFile = new TFile(rootFileName.c_str(), "RECREATE");
0039
0040 theCalibWheel =
0041 pset.getUntrackedParameter<string>("calibWheel", "All");
0042 if (theCalibWheel != "All") {
0043 stringstream linestr;
0044 int selWheel;
0045 linestr << theCalibWheel;
0046 linestr >> selWheel;
0047 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMSPerLayer] chosen wheel " << selWheel;
0048 }
0049
0050
0051 theCalibSector =
0052 pset.getUntrackedParameter<string>("calibSector", "All");
0053 if (theCalibSector != "All") {
0054 stringstream linestr;
0055 int selSector;
0056 linestr << theCalibSector;
0057 linestr >> selSector;
0058 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMSPerLayer] chosen sector " << selSector;
0059 }
0060
0061 vector<string> defaultCell;
0062 defaultCell.push_back("None");
0063 cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell);
0064 for (vector<string>::const_iterator cell = cellsWithHistos.begin(); cell != cellsWithHistos.end(); ++cell) {
0065 if ((*cell) != "None") {
0066 stringstream linestr;
0067 int wheel, sector, station, sl, layer, wire;
0068 linestr << (*cell);
0069 linestr >> wheel >> sector >> station >> sl >> layer >> wire;
0070 wireIdWithHistos.push_back(DTWireId(wheel, station, sector, sl, layer, wire));
0071 }
0072 }
0073
0074 hT0SectorHisto = nullptr;
0075
0076 nevents = 0;
0077 eventsForLayerT0 = pset.getParameter<unsigned int>("eventsForLayerT0");
0078 eventsForWireT0 = pset.getParameter<unsigned int>("eventsForWireT0");
0079 rejectDigiFromPeak = pset.getParameter<unsigned int>("rejectDigiFromPeak");
0080 tpPeakWidth = pset.getParameter<double>("tpPeakWidth");
0081
0082 correctByChamberMean_ = pset.getParameter<bool>("correctByChamberMean");
0083 }
0084
0085
0086 DTT0CalibrationRMS::~DTT0CalibrationRMS() {
0087 if (debug)
0088 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS]Destructor called!";
0089
0090 theFile->Close();
0091 }
0092
0093
0094 void DTT0CalibrationRMS::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0095 if (debug || event.id().event() % 500 == 0)
0096 edm::LogVerbatim("DTCalibration") << "--- [DTT0CalibrationRMS] Analysing Event: #Run: " << event.id().run()
0097 << " #Event: " << event.id().event();
0098 nevents++;
0099
0100
0101 const Handle<DTDigiCollection>& digis = event.getHandle(digiToken);
0102
0103
0104 dtGeom = eventSetup.getHandle(dtGeomToken_);
0105
0106
0107 DTDigiCollection::DigiRangeIterator dtLayerIt;
0108 for (dtLayerIt = digis->begin(); dtLayerIt != digis->end(); ++dtLayerIt) {
0109
0110 const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
0111
0112
0113 const DTLayerId layerId = (*dtLayerIt).first;
0114
0115 if ((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel))
0116 continue;
0117 if ((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector))
0118 continue;
0119
0120
0121 for (DTDigiCollection::const_iterator digi = digiRange.first; digi != digiRange.second; ++digi) {
0122 double t0 = (*digi).countsTDC();
0123
0124
0125 if (nevents < eventsForLayerT0) {
0126
0127 TH1I* hT0LayerHisto = theHistoLayerMap[layerId];
0128
0129 if (hT0LayerHisto == nullptr) {
0130 theFile->cd();
0131 hT0LayerHisto = new TH1I(getHistoName(layerId).c_str(),
0132 "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
0133 200,
0134 t0 - 100,
0135 t0 + 100);
0136 if (debug)
0137 edm::LogVerbatim("DTCalibration") << " New T0 per Layer Histo: " << hT0LayerHisto->GetName();
0138 theHistoLayerMap[layerId] = hT0LayerHisto;
0139 }
0140
0141
0142 theFile->cd();
0143 if (hT0LayerHisto != nullptr) {
0144 hT0LayerHisto->Fill(t0);
0145 }
0146 }
0147
0148
0149 if (nevents > eventsForLayerT0) {
0150
0151 const DTWireId wireId(layerId, (*digi).wire());
0152 if (debug) {
0153 edm::LogVerbatim("DTCalibration")
0154 << " Wire: " << wireId << "\n time (TDC counts): " << (*digi).countsTDC();
0155 }
0156
0157
0158 vector<DTWireId>::iterator it_wire = find(wireIdWithHistos.begin(), wireIdWithHistos.end(), wireId);
0159 if (it_wire != wireIdWithHistos.end()) {
0160 if (theHistoWireMap.find(wireId) == theHistoWireMap.end()) {
0161 theHistoWireMap[wireId] = new TH1I(getHistoName(wireId).c_str(),
0162 "T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",
0163 7000,
0164 0,
0165 7000);
0166 if (debug)
0167 edm::LogVerbatim("DTCalibration") << " New T0 per wire Histo: " << (theHistoWireMap[wireId])->GetName();
0168 }
0169 if (theHistoWireMap_ref.find(wireId) == theHistoWireMap_ref.end()) {
0170 theHistoWireMap_ref[wireId] = new TH1I((getHistoName(wireId) + "_ref").c_str(),
0171 "T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",
0172 7000,
0173 0,
0174 7000);
0175 if (debug)
0176 edm::LogVerbatim("DTCalibration")
0177 << " New T0 per wire Histo: " << (theHistoWireMap_ref[wireId])->GetName();
0178 }
0179
0180 TH1I* hT0WireHisto = theHistoWireMap[wireId];
0181
0182 theFile->cd();
0183 if (hT0WireHisto)
0184 hT0WireHisto->Fill(t0);
0185 }
0186
0187
0188 if (abs(hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin()) - t0) > rejectDigiFromPeak) {
0189 if (debug)
0190 edm::LogVerbatim("DTCalibration") << "digi skipped because t0 per sector "
0191 << hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin());
0192 continue;
0193 }
0194
0195
0196 if (nevents < (eventsForLayerT0 + eventsForWireT0)) {
0197
0198 if (it_wire != wireIdWithHistos.end()) {
0199 TH1I* hT0WireHisto_ref = theHistoWireMap_ref[wireId];
0200 theFile->cd();
0201 if (hT0WireHisto_ref)
0202 hT0WireHisto_ref->Fill(t0);
0203 }
0204 if (!nDigiPerWire_ref[wireId]) {
0205 mK_ref[wireId] = 0;
0206 }
0207 nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1;
0208 mK_ref[wireId] = mK_ref[wireId] + (t0 - mK_ref[wireId]) / nDigiPerWire_ref[wireId];
0209 }
0210
0211 else if (nevents > (eventsForLayerT0 + eventsForWireT0)) {
0212 if (abs(t0 - mK_ref[wireId]) > tpPeakWidth)
0213 continue;
0214 if (!nDigiPerWire[wireId]) {
0215 theAbsoluteT0PerWire[wireId] = 0;
0216 qK[wireId] = 0;
0217 mK[wireId] = 0;
0218 }
0219 nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1;
0220 theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0;
0221
0222 qK[wireId] =
0223 qK[wireId] + ((nDigiPerWire[wireId] - 1) * (t0 - mK[wireId]) * (t0 - mK[wireId]) / nDigiPerWire[wireId]);
0224 mK[wireId] = mK[wireId] + (t0 - mK[wireId]) / nDigiPerWire[wireId];
0225 }
0226 }
0227 }
0228 }
0229
0230
0231 if (nevents == eventsForLayerT0) {
0232 for (map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin(); lHisto != theHistoLayerMap.end();
0233 ++lHisto) {
0234 if (debug)
0235 edm::LogVerbatim("DTCalibration") << "Reading histogram " << (*lHisto).second->GetName() << " with mean "
0236 << (*lHisto).second->GetMean() << " and RMS " << (*lHisto).second->GetRMS();
0237
0238
0239 if ((*lHisto).second->GetRMS() < 5.0) {
0240 if (hT0SectorHisto == nullptr) {
0241 hT0SectorHisto = new TH1D("hT0AllLayerOfSector",
0242 "T0 from pulses per layer in sector",
0243
0244 700,
0245 0,
0246 7000);
0247 }
0248 if (debug)
0249 edm::LogVerbatim("DTCalibration") << " accepted";
0250 hT0SectorHisto->Fill((*lHisto).second->GetMean());
0251 }
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267 theT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetMean();
0268 theSigmaT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetRMS();
0269 }
0270 if (!hT0SectorHisto) {
0271 edm::LogVerbatim("DTCalibration")
0272 << "[DTT0CalibrationRMS]: All the t0 per layer are still uncorrect: trying with greater number of events";
0273 eventsForLayerT0 = eventsForLayerT0 * 2;
0274 return;
0275 }
0276 if (debug)
0277 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] t0 reference for this sector "
0278 << hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin());
0279 }
0280 }
0281
0282 void DTT0CalibrationRMS::endJob() {
0283 DTT0* t0sAbsolute = new DTT0();
0284 DTT0* t0sRelative = new DTT0();
0285 DTT0* t0sWRTChamber = new DTT0();
0286
0287
0288 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMSPerLayer]Writing histos to file!";
0289
0290 theFile->cd();
0291 theFile->WriteTObject(hT0SectorHisto);
0292
0293 for (map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap.begin(); wHisto != theHistoWireMap.end();
0294 ++wHisto) {
0295 (*wHisto).second->Write();
0296 }
0297 for (map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap_ref.begin(); wHisto != theHistoWireMap_ref.end();
0298 ++wHisto) {
0299 (*wHisto).second->Write();
0300 }
0301 for (map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin(); lHisto != theHistoLayerMap.end();
0302 ++lHisto) {
0303 (*lHisto).second->Write();
0304 }
0305
0306
0307 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] Compute and store t0 and sigma per wire";
0308
0309 for (map<DTWireId, double>::const_iterator wiret0 = theAbsoluteT0PerWire.begin();
0310 wiret0 != theAbsoluteT0PerWire.end();
0311 ++wiret0) {
0312 if (nDigiPerWire[(*wiret0).first]) {
0313 double t0 = (*wiret0).second / nDigiPerWire[(*wiret0).first];
0314
0315 theRelativeT0PerWire[(*wiret0).first] = t0 - hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin());
0316
0317
0318 theSigmaT0PerWire[(*wiret0).first] = sqrt(qK[(*wiret0).first] / nDigiPerWire[(*wiret0).first]);
0319
0320 edm::LogVerbatim("DTCalibration") << "Wire " << (*wiret0).first << " has t0 " << t0 << "(absolute) "
0321 << theRelativeT0PerWire[(*wiret0).first] << "(relative)"
0322 << " sigma " << theSigmaT0PerWire[(*wiret0).first];
0323
0324 t0sAbsolute->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first], DTTimeUnits::counts);
0325 } else {
0326 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] ERROR: no digis in wire " << (*wiret0).first;
0327 abort();
0328 }
0329 }
0330
0331 if (correctByChamberMean_) {
0332
0333
0334 const vector<const DTSuperLayer*> superLayers = dtGeom->superLayers();
0335
0336 for (vector<const DTSuperLayer*>::const_iterator sl = superLayers.begin(); sl != superLayers.end(); sl++) {
0337
0338 double oddLayersMean = 0;
0339 double evenLayersMean = 0;
0340 double oddLayersDen = 0;
0341 double evenLayersDen = 0;
0342 for (map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
0343 wiret0 != theRelativeT0PerWire.end();
0344 ++wiret0) {
0345 if ((*wiret0).first.layerId().superlayerId() == (*sl)->id()) {
0346 if (debug)
0347 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] Superlayer " << (*sl)->id() << "layer "
0348 << (*wiret0).first.layerId().layer() << " with " << (*wiret0).second;
0349 if (((*wiret0).first.layerId().layer()) % 2) {
0350 oddLayersMean = oddLayersMean + (*wiret0).second;
0351 oddLayersDen++;
0352 } else {
0353 evenLayersMean = evenLayersMean + (*wiret0).second;
0354 evenLayersDen++;
0355 }
0356 }
0357 }
0358 oddLayersMean = oddLayersMean / oddLayersDen;
0359 evenLayersMean = evenLayersMean / evenLayersDen;
0360
0361 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] Relative T0 mean for odd layers " << oddLayersMean
0362 << " even layers" << evenLayersMean;
0363
0364
0365 double oddLayersSigma = 0;
0366 double evenLayersSigma = 0;
0367 for (map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
0368 wiret0 != theRelativeT0PerWire.end();
0369 ++wiret0) {
0370 if ((*wiret0).first.layerId().superlayerId() == (*sl)->id()) {
0371 if (((*wiret0).first.layerId().layer()) % 2) {
0372 oddLayersSigma = oddLayersSigma + ((*wiret0).second - oddLayersMean) * ((*wiret0).second - oddLayersMean);
0373 } else {
0374 evenLayersSigma =
0375 evenLayersSigma + ((*wiret0).second - evenLayersMean) * ((*wiret0).second - evenLayersMean);
0376 }
0377 }
0378 }
0379 oddLayersSigma = oddLayersSigma / oddLayersDen;
0380 evenLayersSigma = evenLayersSigma / evenLayersDen;
0381 oddLayersSigma = sqrt(oddLayersSigma);
0382 evenLayersSigma = sqrt(evenLayersSigma);
0383
0384
0385 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] Relative T0 sigma for odd layers " << oddLayersSigma
0386 << " even layers" << evenLayersSigma;
0387
0388
0389 double oddLayersFinalMean = 0;
0390 double evenLayersFinalMean = 0;
0391 for (map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
0392 wiret0 != theRelativeT0PerWire.end();
0393 ++wiret0) {
0394 if ((*wiret0).first.layerId().superlayerId() == (*sl)->id()) {
0395 if (((*wiret0).first.layerId().layer()) % 2) {
0396 if (abs((*wiret0).second - oddLayersMean) < (2 * oddLayersSigma))
0397 oddLayersFinalMean = oddLayersFinalMean + (*wiret0).second;
0398 } else {
0399 if (abs((*wiret0).second - evenLayersMean) < (2 * evenLayersSigma))
0400 evenLayersFinalMean = evenLayersFinalMean + (*wiret0).second;
0401 }
0402 }
0403 }
0404 oddLayersFinalMean = oddLayersFinalMean / oddLayersDen;
0405 evenLayersFinalMean = evenLayersFinalMean / evenLayersDen;
0406
0407 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] Final relative T0 mean for odd layers "
0408 << oddLayersFinalMean << " even layers" << evenLayersFinalMean;
0409
0410 for (map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
0411 wiret0 != theRelativeT0PerWire.end();
0412 ++wiret0) {
0413 if ((*wiret0).first.layerId().superlayerId() == (*sl)->id()) {
0414 double t0 = -999;
0415 if (((*wiret0).first.layerId().layer()) % 2)
0416 t0 = (*wiret0).second + (evenLayersFinalMean - oddLayersFinalMean);
0417 else
0418 t0 = (*wiret0).second;
0419
0420 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS] Wire " << (*wiret0).first << " has t0 "
0421 << (*wiret0).second << " (relative, after even-odd layer corrections) "
0422 << " sigma " << theSigmaT0PerWire[(*wiret0).first];
0423
0424
0425 t0sRelative->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first], DTTimeUnits::counts);
0426 }
0427 }
0428 }
0429
0430
0431 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS]Computing relative t0 wrt to chamber average";
0432
0433 map<DTChamberId, double> sumT0ByChamber;
0434 map<DTChamberId, int> countT0ByChamber;
0435 for (DTT0::const_iterator tzero = t0sRelative->begin(); tzero != t0sRelative->end(); ++tzero) {
0436 int channelId = tzero->channelId;
0437 if (channelId == 0)
0438 continue;
0439 DTWireId wireId(channelId);
0440 DTChamberId chamberId(wireId.chamberId());
0441
0442
0443 float t0mean_f;
0444 float t0rms_f;
0445 t0sRelative->get(wireId, t0mean_f, t0rms_f, DTTimeUnits::counts);
0446 sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + t0mean_f;
0447
0448 countT0ByChamber[chamberId]++;
0449 }
0450
0451
0452 for (DTT0::const_iterator tzero = t0sRelative->begin(); tzero != t0sRelative->end(); ++tzero) {
0453 int channelId = tzero->channelId;
0454 if (channelId == 0)
0455 continue;
0456 DTWireId wireId(channelId);
0457 DTChamberId chamberId(wireId.chamberId());
0458
0459
0460
0461 float t0mean_f;
0462 float t0rms_f;
0463 t0sRelative->get(wireId, t0mean_f, t0rms_f, DTTimeUnits::counts);
0464 double t0mean = t0mean_f - (sumT0ByChamber[chamberId] / countT0ByChamber[chamberId]);
0465 double t0rms = t0rms_f;
0466
0467 t0sWRTChamber->set(wireId, t0mean, t0rms, DTTimeUnits::counts);
0468 edm::LogVerbatim("DTCalibration") << "Changing t0 of wire " << wireId << " from " << t0mean_f << " to " << t0mean;
0469 }
0470 }
0471
0472
0473 if (debug)
0474 edm::LogVerbatim("DTCalibration") << "[DTT0CalibrationRMS]Writing values in DB!";
0475
0476 string t0Record = "DTT0Rcd";
0477
0478 if (correctByChamberMean_)
0479 DTCalibDBUtils::writeToDB(t0Record, t0sWRTChamber);
0480 else
0481 DTCalibDBUtils::writeToDB(t0Record, t0sAbsolute);
0482 }
0483
0484 string DTT0CalibrationRMS::getHistoName(const DTWireId& wId) const {
0485 string histoName = "Ch_" + std::to_string(wId.wheel()) + "_" + std::to_string(wId.station()) + "_" +
0486 std::to_string(wId.sector()) + "_SL" + std::to_string(wId.superlayer()) + "_L" +
0487 std::to_string(wId.layer()) + "_W" + std::to_string(wId.wire()) + "_hT0Histo";
0488 return histoName;
0489 }
0490
0491 string DTT0CalibrationRMS::getHistoName(const DTLayerId& lId) const {
0492 string histoName = "Ch_" + std::to_string(lId.wheel()) + "_" + std::to_string(lId.station()) + "_" +
0493 std::to_string(lId.sector()) + "_SL" + std::to_string(lId.superlayer()) + "_L" +
0494 std::to_string(lId.layer()) + "_hT0Histo";
0495 return histoName;
0496 }