DTCalibration/plugins/DTTTrigCalibration.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Cerminara - INFN Torino
0005  */
0006 #include "CalibMuon/DTCalibration/plugins/DTTTrigCalibration.h"
0007 #include "CalibMuon/DTCalibration/interface/DTTimeBoxFitter.h"
0008 #include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
0009 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0010
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016
0017 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0018
0019 #include "CondFormats/DTObjects/interface/DTTtrig.h"
0020
0021 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0022
0023 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0024
0025 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
0026 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0027
0028 #include "TFile.h"
0029 #include "TH1F.h"
0030 #include "TGraph.h"
0031
0032 class DTLayer;
0033
0034 using namespace std;
0035 using namespace edm;
0036 // using namespace cond;
0037
0038 // Constructor
0039 DTTTrigCalibration::DTTTrigCalibration(const edm::ParameterSet& pset) {
0040   // Get the debug parameter for verbose output
0041   debug = pset.getUntrackedParameter<bool>("debug");
0042
0043   // Get the label to retrieve digis from the event
0044   std::string digiLabel = pset.getUntrackedParameter<string>("digiLabel");
0045   digiToken = consumes<DTDigiCollection>(edm::InputTag(digiLabel));
0046
0047   // Switch on/off the DB writing
0048   findTMeanAndSigma = pset.getUntrackedParameter<bool>("fitAndWrite", false);
0049
0050   // The TDC time-window (ns)
0051   maxTDCCounts = 5000 * pset.getUntrackedParameter<int>("tdcRescale", 1);
0052   //The maximum number of digis per layer
0053   maxDigiPerLayer = pset.getUntrackedParameter<int>("maxDigiPerLayer", 10);
0054
0055   // The root file which will contain the histos
0056   string rootFileName = pset.getUntrackedParameter<string>("rootFileName");
0057   theFile = new TFile(rootFileName.c_str(), "RECREATE");
0058   theFile->cd();
0059   theFitter = std::make_unique<DTTimeBoxFitter>();
0060   if (debug)
0061     theFitter->setVerbosity(1);
0062
0063   double sigmaFit = pset.getUntrackedParameter<double>("sigmaTTrigFit", 10.);
0064   theFitter->setFitSigma(sigmaFit);
0065
0066   doSubtractT0 = pset.getUntrackedParameter<bool>("doSubtractT0", "false");
0067   // Get the synchronizer
0068   if (doSubtractT0) {
0069     theSync = DTTTrigSyncFactory::get()->create(pset.getUntrackedParameter<string>("tTrigMode"),
0070                                                 pset.getUntrackedParameter<ParameterSet>("tTrigModeConfig"),
0071                                                 consumesCollector());
0072   }
0073
0074   checkNoisyChannels = pset.getUntrackedParameter<bool>("checkNoisyChannels", "false");
0075
0076   // the kfactor to be uploaded in the ttrig DB
0077   kFactor = pset.getUntrackedParameter<double>("kFactor", -0.7);
0078
0079   if (debug)
0080     cout << "[DTTTrigCalibration]Constructor called!" << endl;
0081
0082   if (checkNoisyChannels) {
0083     theStatusMapToken = esConsumes();
0084   }
0085 }
0086
0087 // Destructor
0088 DTTTrigCalibration::~DTTTrigCalibration() {
0089   if (debug)
0090     cout << "[DTTTrigCalibration]Destructor called!" << endl;
0091
0092   //   // Delete all histos
0093   //   for(map<DTSuperLayerId, TH1F*>::const_iterator slHisto = theHistoMap.begin();
0094   //       slHisto != theHistoMap.end();
0095   //       slHisto++) {
0096   //     delete (*slHisto).second;
0097   //   }
0098
0099   theFile->Close();
0100 }
0101
0102 /// Perform the real analysis
0103 void DTTTrigCalibration::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0104   if (debug)
0105     cout << "[DTTTrigCalibration] #Event: " << event.id().event() << endl;
0106
0107   // Get the digis from the event
0108   const edm::Handle<DTDigiCollection>& digis = event.getHandle(digiToken);
0109
0110   ESHandle<DTStatusFlag> statusMap;
0111   if (checkNoisyChannels) {
0112     // Get the map of noisy channels
0113     statusMap = eventSetup.getHandle(theStatusMapToken);
0114   }
0115
0116   if (doSubtractT0)
0117     theSync->setES(eventSetup);
0118
0119   //The chambers too noisy in this event
0120   vector<DTChamberId> badChambers;
0121
0122   // Iterate through all digi collections ordered by LayerId
0123   DTDigiCollection::DigiRangeIterator dtLayerIt;
0124   for (dtLayerIt = digis->begin(); dtLayerIt != digis->end(); ++dtLayerIt) {
0125     // Get the iterators over the digis associated with this LayerId
0126     const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
0127
0128     const DTLayerId layerId = (*dtLayerIt).first;
0129     const DTSuperLayerId slId = layerId.superlayerId();
0130     const DTChamberId chId = slId.chamberId();
0131     bool badChamber = false;
0132
0133     if (debug)
0134       cout << "----------- Layer " << layerId << " -------------" << endl;
0135
0136     //Check if the layer is inside a noisy chamber
0137     for (vector<DTChamberId>::const_iterator chamber = badChambers.begin(); chamber != badChambers.end(); ++chamber) {
0138       if ((*chamber) == chId) {
0139         badChamber = true;
0140         break;
0141       }
0142     }
0143     if (badChamber)
0144       continue;
0145
0146     //Check if the layer has too many digis
0147     if ((digiRange.second - digiRange.first) > maxDigiPerLayer) {
0148       if (debug)
0149         cout << "Layer " << layerId << "has too many digis (" << (digiRange.second - digiRange.first) << ")" << endl;
0150       badChambers.push_back(chId);
0151       continue;
0152     }
0153
0154     // Get the histo from the map
0155     TH1F* hTBox = theHistoMap[slId];
0156     if (hTBox == nullptr) {
0157       // Book the histogram
0158       theFile->cd();
0159       hTBox =
0160           new TH1F(getTBoxName(slId).c_str(), "Time box (ns)", int(0.25 * 32.0 * maxTDCCounts / 25.0), 0, maxTDCCounts);
0161       if (debug)
0162         cout << "  New Time Box: " << hTBox->GetName() << endl;
0163       theHistoMap[slId] = hTBox;
0164     }
0165     TH1F* hO = theOccupancyMap[layerId];
0166     if (hO == nullptr) {
0167       // Book the histogram
0168       theFile->cd();
0169       hO = new TH1F(getOccupancyName(layerId).c_str(), "Occupancy", 100, 0, 100);
0170       if (debug)
0171         cout << "  New Time Box: " << hO->GetName() << endl;
0172       theOccupancyMap[layerId] = hO;
0173     }
0174
0175     // Loop over all digis in the given range
0176     for (DTDigiCollection::const_iterator digi = digiRange.first; digi != digiRange.second; ++digi) {
0177       const DTWireId wireId(layerId, (*digi).wire());
0178
0179       // Check for noisy channels and skip them
0180       if (checkNoisyChannels) {
0181         bool isNoisy = false;
0182         bool isFEMasked = false;
0183         bool isTDCMasked = false;
0184         bool isTrigMask = false;
0185         bool isDead = false;
0186         bool isNohv = false;
0187         statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
0188         if (isNoisy) {
0189           if (debug)
0190             cout << "Wire: " << wireId << " is noisy, skipping!" << endl;
0191           continue;
0192         }
0193       }
0194       theFile->cd();
0195       double offset = 0;
0196       if (doSubtractT0) {
0197         const DTLayer* layer = nullptr;  //fake
0198         const GlobalPoint glPt;          //fake
0199         offset = theSync->offset(layer, wireId, glPt);
0200       }
0201       hTBox->Fill((*digi).time() - offset);
0202       if (debug) {
0203         cout << "   Filling Time Box:   " << hTBox->GetName() << endl;
0204         cout << "           offset (ns): " << offset << endl;
0205         cout << "           time(ns):   " << (*digi).time() - offset << endl;
0206       }
0207       hO->Fill((*digi).wire());
0208     }
0209   }
0210 }
0211
0212 void DTTTrigCalibration::endJob() {
0213   if (debug)
0214     cout << "[DTTTrigCalibration]Writing histos to file!" << endl;
0215
0216   // Write all time boxes to file
0217   theFile->cd();
0218   for (map<DTSuperLayerId, TH1F*>::const_iterator slHisto = theHistoMap.begin(); slHisto != theHistoMap.end();
0219        ++slHisto) {
0220     (*slHisto).second->Write();
0221   }
0222   for (map<DTLayerId, TH1F*>::const_iterator slHisto = theOccupancyMap.begin(); slHisto != theOccupancyMap.end();
0223        ++slHisto) {
0224     (*slHisto).second->Write();
0225   }
0226
0227   if (findTMeanAndSigma) {
0228     // Create the object to be written to DB
0229     DTTtrig* tTrig = new DTTtrig();
0230
0231     // Loop over the map, fit the histos and write the resulting values to the DB
0232     for (map<DTSuperLayerId, TH1F*>::const_iterator slHisto = theHistoMap.begin(); slHisto != theHistoMap.end();
0233          ++slHisto) {
0234       pair<double, double> meanAndSigma = theFitter->fitTimeBox((*slHisto).second);
0235       tTrig->set((*slHisto).first, meanAndSigma.first, meanAndSigma.second, kFactor, DTTimeUnits::ns);
0236
0237       if (debug) {
0238         cout << " SL: " << (*slHisto).first << " mean = " << meanAndSigma.first << " sigma = " << meanAndSigma.second
0239              << endl;
0240       }
0241     }
0242
0243     // Print the ttrig map
0244     dumpTTrigMap(tTrig);
0245
0246     // Plot the tTrig
0247     plotTTrig(tTrig);
0248
0249     if (debug)
0250       cout << "[DTTTrigCalibration]Writing ttrig object to DB!" << endl;
0251
0252     // FIXME: to be read from cfg?
0253     string tTrigRecord = "DTTtrigRcd";
0254
0255     // Write the object to DB
0256     DTCalibDBUtils::writeToDB(tTrigRecord, tTrig);
0257   }
0258 }
0259
0260 string DTTTrigCalibration::getTBoxName(const DTSuperLayerId& slId) const {
0261   string histoName;
0262   stringstream theStream;
0263   theStream << "Ch_" << slId.wheel() << "_" << slId.station() << "_" << slId.sector() << "_SL" << slId.superlayer()
0264             << "_hTimeBox";
0265   theStream >> histoName;
0266   return histoName;
0267 }
0268
0269 string DTTTrigCalibration::getOccupancyName(const DTLayerId& slId) const {
0270   string histoName;
0271   stringstream theStream;
0272   theStream << "Ch_" << slId.wheel() << "_" << slId.station() << "_" << slId.sector() << "_SL" << slId.superlayer()
0273             << "_L" << slId.layer() << "_Occupancy";
0274   theStream >> histoName;
0275   return histoName;
0276 }
0277
0278 void DTTTrigCalibration::dumpTTrigMap(const DTTtrig* tTrig) const {
0279   static const double convToNs = 25. / 32.;
0280   for (DTTtrig::const_iterator ttrig = tTrig->begin(); ttrig != tTrig->end(); ++ttrig) {
0281     cout << "Wh: " << (*ttrig).first.wheelId << " St: " << (*ttrig).first.stationId
0282          << " Sc: " << (*ttrig).first.sectorId << " Sl: " << (*ttrig).first.slId
0283          << " TTrig mean (ns): " << (*ttrig).second.tTrig * convToNs
0284          << " TTrig sigma (ns): " << (*ttrig).second.tTrms * convToNs << endl;
0285   }
0286 }
0287
0288 void DTTTrigCalibration::plotTTrig(const DTTtrig* tTrig) const {
0289   TH1F* tTrig_YB1_Se10 = new TH1F("tTrig_YB1_Se10", "tTrig YB1_Se10", 15, 1, 16);
0290   TH1F* tTrig_YB2_Se10 = new TH1F("tTrig_YB2_Se10", "tTrig YB2_Se10", 15, 1, 16);
0291   TH1F* tTrig_YB2_Se11 = new TH1F("tTrig_YB2_Se11", "tTrig YB2_Se11", 12, 1, 13);
0292
0293   static const double convToNs = 25. / 32.;
0294   for (DTTtrig::const_iterator ttrig = tTrig->begin(); ttrig != tTrig->end(); ++ttrig) {
0295     // avoid to have wired numbers in the plot
0296     float tTrigValue = 0;
0297     float tTrmsValue = 0;
0298     if ((*ttrig).second.tTrig * convToNs > 0 && (*ttrig).second.tTrig * convToNs < 32000) {
0299       tTrigValue = (*ttrig).second.tTrig * convToNs;
0300       tTrmsValue = (*ttrig).second.tTrms * convToNs;
0301     }
0302
0303     int binx;
0304     string binLabel;
0305     stringstream binLabelStream;
0306     if ((*ttrig).first.sectorId != 14) {
0307       binx = ((*ttrig).first.stationId - 1) * 3 + (*ttrig).first.slId;
0308       binLabelStream << "MB" << (*ttrig).first.stationId << "_SL" << (*ttrig).first.slId;
0309     } else {
0310       binx = 12 + (*ttrig).first.slId;
0311       binLabelStream << "MB14_SL" << (*ttrig).first.slId;
0312     }
0313     binLabelStream >> binLabel;
0314
0315     if ((*ttrig).first.wheelId == 2) {
0316       if ((*ttrig).first.sectorId == 10 || (*ttrig).first.sectorId == 14) {
0317         tTrig_YB2_Se10->Fill(binx, tTrigValue);
0318         tTrig_YB2_Se10->SetBinError(binx, tTrmsValue);
0319         tTrig_YB2_Se10->GetXaxis()->SetBinLabel(binx, binLabel.c_str());
0320         tTrig_YB2_Se10->GetYaxis()->SetTitle("ns");
0321       } else {
0322         tTrig_YB2_Se11->Fill(binx, tTrigValue);
0323         tTrig_YB2_Se11->SetBinError(binx, tTrmsValue);
0324         tTrig_YB2_Se11->GetXaxis()->SetBinLabel(binx, binLabel.c_str());
0325         tTrig_YB2_Se11->GetYaxis()->SetTitle("ns");
0326       }
0327     } else {
0328       tTrig_YB1_Se10->Fill(binx, tTrigValue);
0329       tTrig_YB1_Se10->SetBinError(binx, tTrmsValue);
0330       tTrig_YB1_Se10->GetXaxis()->SetBinLabel(binx, binLabel.c_str());
0331       tTrig_YB1_Se10->GetYaxis()->SetTitle("ns");
0332     }
0333   }
0334
0335   tTrig_YB1_Se10->Write();
0336   tTrig_YB2_Se10->Write();
0337   tTrig_YB2_Se11->Write();
0338 }