Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:34:56 

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  $Date: 2012/05/11 17:17:17 $
0005  *  $Revision: 1.6 $
0006  *  \author S. Bolognesi - INFN Torino
0007  *  06/08/2008 Mofified by Antonio.Vilela.Pereira@cern.ch
0008  */
0009 
0010 #include "CalibMuon/DTCalibration/plugins/DTT0Calibration.h"
0011 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0012 
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 
0016 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0017 #include "Geometry/Records/interface/MuonNumberingRecord.h"
0018 
0019 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0020 #include "CondFormats/DTObjects/interface/DTT0.h"
0021 
0022 #include <CondFormats/DTObjects/interface/DTTtrig.h>
0023 #include <CondFormats/DataRecord/interface/DTTtrigRcd.h>
0024 
0025 #include "TKey.h"
0026 #include "TF1.h"
0027 
0028 #include <cassert>
0029 
0030 using namespace std;
0031 using namespace edm;
0032 
0033 // Constructor
0034 DTT0Calibration::DTT0Calibration(const edm::ParameterSet& pset)
0035     : debug(pset.getUntrackedParameter<bool>("debug")),
0036       digiToken(consumes<DTDigiCollection>(pset.getUntrackedParameter<string>("digiLabel"))),
0037       theFile(pset.getUntrackedParameter<string>("rootFileName", "DTT0PerLayer.root").c_str(), "RECREATE"),
0038       nevents(0),
0039       eventsForLayerT0(pset.getParameter<unsigned int>("eventsForLayerT0")),
0040       eventsForWireT0(pset.getParameter<unsigned int>("eventsForWireT0")),
0041       tpPeakWidth(pset.getParameter<double>("tpPeakWidth")),
0042       tpPeakWidthPerLayer(pset.getParameter<double>("tpPeakWidthPerLayer")),
0043       rejectDigiFromPeak(pset.getParameter<unsigned int>("rejectDigiFromPeak")),
0044       hLayerPeaks("hLayerPeaks", "", 3000, 0, 3000),
0045       spectrum(20),
0046       dtGeomToken_(esConsumes()) {
0047   // Get the debug parameter for verbose output
0048   if (debug)
0049     cout << "[DTT0Calibration]Constructor called!" << endl;
0050 
0051   theCalibWheel =
0052       pset.getUntrackedParameter<string>("calibWheel", "All");  //FIXME amke a vector of integer instead of a string
0053   if (theCalibWheel != "All") {
0054     stringstream linestr;
0055     int selWheel;
0056     linestr << theCalibWheel;
0057     linestr >> selWheel;
0058     cout << "[DTT0CalibrationPerLayer] chosen wheel " << selWheel << endl;
0059   }
0060 
0061   // Sector/s to calibrate
0062   theCalibSector =
0063       pset.getUntrackedParameter<string>("calibSector", "All");  //FIXME amke a vector of integer instead of a string
0064   if (theCalibSector != "All") {
0065     stringstream linestr;
0066     int selSector;
0067     linestr << theCalibSector;
0068     linestr >> selSector;
0069     cout << "[DTT0CalibrationPerLayer] chosen sector " << selSector << endl;
0070   }
0071 
0072   vector<string> defaultCell;
0073   const auto& cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell);
0074   for (const auto& cell : cellsWithHistos) {
0075     stringstream linestr;
0076     int wheel, sector, station, sl, layer, wire;
0077     linestr << cell;
0078     linestr >> wheel >> sector >> station >> sl >> layer >> wire;
0079     wireIdWithHistos.push_back(DTWireId(wheel, station, sector, sl, layer, wire));
0080   }
0081 }
0082 
0083 // Destructor
0084 DTT0Calibration::~DTT0Calibration() {
0085   if (debug)
0086     cout << "[DTT0Calibration]Destructor called!" << endl;
0087 
0088   theFile.Close();
0089 }
0090 
0091 /// Perform the real analysis
0092 void DTT0Calibration::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0093   nevents++;
0094 
0095   // Get the digis from the event
0096   Handle<DTDigiCollection> digis;
0097   event.getByToken(digiToken, digis);
0098 
0099   // Get the DT Geometry
0100   if (nevents == 1)
0101     dtGeom = eventSetup.getHandle(dtGeomToken_);
0102 
0103   // Iterate through all digi collections ordered by LayerId
0104   for (const auto& digis_per_layer : *digis) {
0105     //std::cout << __LINE__ << std::endl;
0106     // Get the iterators over the digis associated with this LayerId
0107     const DTDigiCollection::Range& digiRange = digis_per_layer.second;
0108 
0109     // Get the layerId
0110     const DTLayerId layerId = digis_per_layer.first;
0111     //const DTChamberId chamberId = layerId.superlayerId().chamberId();
0112 
0113     if ((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel))
0114       continue;
0115     if ((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector))
0116       continue;
0117 
0118     // Loop over all digis in the given layer
0119     for (DTDigiCollection::const_iterator digi = digiRange.first; digi != digiRange.second; ++digi) {
0120       const double t0 = digi->countsTDC();
0121       const DTWireId wireIdtmp(layerId, (*digi).wire());
0122 
0123       // Use first bunch of events to fill t0 per layer
0124       if (nevents <= eventsForLayerT0) {
0125         // If it doesn't exist, book it
0126         if (not theHistoLayerMap.count(layerId)) {
0127           theHistoLayerMap[layerId] = TH1I(getHistoName(layerId).c_str(),
0128                                            "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
0129                                            3000,
0130                                            0,
0131                                            3000);
0132           if (debug)
0133             cout << "  New T0 per Layer Histo: " << theHistoLayerMap[layerId].GetName() << endl;
0134         }
0135         theHistoLayerMap[layerId].Fill(t0);
0136       }
0137 
0138       // Use all the remaining events to compute t0 per wire
0139       if (nevents > eventsForLayerT0) {
0140         // Get the wireId
0141         const DTWireId wireId(layerId, (*digi).wire());
0142         if (debug) {
0143           cout << "   Wire: " << wireId << endl << "       time (TDC counts): " << (*digi).countsTDC() << endl;
0144         }
0145 
0146         //Fill the histos per wire for the chosen cells
0147         if (std::find(layerIdWithWireHistos.begin(), layerIdWithWireHistos.end(), layerId) !=
0148                 layerIdWithWireHistos.end() or
0149             std::find(wireIdWithHistos.begin(), wireIdWithHistos.end(), wireId) != wireIdWithHistos.end()) {
0150           //If it doesn't exist, book it
0151           if (theHistoWireMap.count(wireId) == 0) {
0152             theHistoWireMap[wireId] = TH1I(getHistoName(wireId).c_str(),
0153                                            "T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",
0154                                            7000,
0155                                            0,
0156                                            7000);
0157             if (debug)
0158               cout << "  New T0 per wire Histo: " << theHistoWireMap[wireId].GetName() << endl;
0159           }
0160           theHistoWireMap[wireId].Fill(t0);
0161         }
0162 
0163         //Select per layer
0164         if (fabs(theTPPeakMap[layerId] - t0) > rejectDigiFromPeak) {
0165           if (debug)
0166             cout << "digi skipped because t0 too far from peak " << theTPPeakMap[layerId] << endl;
0167           continue;
0168         }
0169 
0170         //Use second bunch of events to compute a t0 reference per wire
0171         if (nevents <= (eventsForLayerT0 + eventsForWireT0)) {
0172           if (!nDigiPerWire_ref[wireId]) {
0173             mK_ref[wireId] = 0;
0174           }
0175           nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1;
0176           mK_ref[wireId] = mK_ref[wireId] + (t0 - mK_ref[wireId]) / nDigiPerWire_ref[wireId];
0177         }
0178         //Use last all the remaining events to compute the mean and sigma t0 per wire
0179         else if (nevents > (eventsForLayerT0 + eventsForWireT0)) {
0180           if (abs(t0 - mK_ref[wireId]) > tpPeakWidth)
0181             continue;
0182           if (!nDigiPerWire[wireId]) {
0183             theAbsoluteT0PerWire[wireId] = 0;
0184             qK[wireId] = 0;
0185             mK[wireId] = 0;
0186           }
0187           nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1;
0188           theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0;
0189           qK[wireId] =
0190               qK[wireId] + ((nDigiPerWire[wireId] - 1) * (t0 - mK[wireId]) * (t0 - mK[wireId]) / nDigiPerWire[wireId]);
0191           mK[wireId] = mK[wireId] + (t0 - mK[wireId]) / nDigiPerWire[wireId];
0192         }
0193       }  //end if(nevents>1000)
0194     }  //end loop on digi
0195   }  //end loop on layer
0196 
0197   //Use the t0 per layer histos to have an indication about the t0 position
0198   if (nevents == eventsForLayerT0) {
0199     for (const auto& lHisto : theHistoLayerMap) {
0200       const auto& layerId = lHisto.first;
0201       const auto& hist = lHisto.second;
0202       if (debug)
0203         cout << "Reading histogram " << hist.GetName() << " with mean " << hist.GetMean() << " and RMS "
0204              << hist.GetRMS() << endl;
0205 
0206       //Find peaks
0207       int npeaks = spectrum.Search(&hist, (tpPeakWidthPerLayer / 2.), "", 0.3);
0208 
0209       double* peaks = spectrum.GetPositionX();
0210       //Put in a std::vector<float>
0211       vector<double> peakMeans(peaks, peaks + npeaks);
0212       //Sort the peaks in ascending order
0213       sort(peakMeans.begin(), peakMeans.end());
0214 
0215       if (peakMeans.empty()) {
0216         theTPPeakMap[layerId] = hist.GetMaximumBin();
0217         std::cout << "No peaks found by peakfinder in layer " << layerId << ". Taking maximum bin at "
0218                   << theTPPeakMap[layerId] << ". Please check!" << std::endl;
0219         layerIdWithWireHistos.push_back(layerId);
0220       } else if (fabs(hist.GetXaxis()->FindBin(peakMeans.front()) - hist.GetXaxis()->FindBin(peakMeans.back())) <
0221                  rejectDigiFromPeak) {
0222         theTPPeakMap[layerId] = peakMeans[peakMeans.size() / 2];
0223       } else {
0224         bool peak_set = false;
0225         for (const auto& peak : peakMeans) {
0226           // Skip if at low edge
0227           if (peak - tpPeakWidthPerLayer <= 0)
0228             continue;
0229           // Get integral of peak
0230           double sum = 0;
0231           for (int ibin = peak - tpPeakWidthPerLayer; ibin < peak + tpPeakWidthPerLayer; ibin++) {
0232             sum += hist.GetBinContent(ibin);
0233           }
0234           // Skip if peak too small
0235           if (sum < hist.GetMaximum() / 2)
0236             continue;
0237 
0238           // Passed all cuts
0239           theTPPeakMap[layerId] = peak;
0240           peak_set = true;
0241           break;
0242         }
0243         if (peak_set) {
0244           std::cout << "Peaks to far away from each other in layer " << layerId
0245                     << ". Maybe cross talk? Taking first good peak at " << theTPPeakMap[layerId] << ". Please check!"
0246                     << std::endl;
0247           layerIdWithWireHistos.push_back(layerId);
0248         } else {
0249           theTPPeakMap[layerId] = hist.GetMaximumBin();
0250           std::cout << "Peaks to far away from each other in layer " << layerId
0251                     << " and no good peak found. Taking maximum bin at " << theTPPeakMap[layerId] << ". Please check!"
0252                     << std::endl;
0253           layerIdWithWireHistos.push_back(layerId);
0254         }
0255       }
0256       if (peakMeans.size() > 5) {
0257         std::cout << "Found more than 5 peaks in layer " << layerId << ". Please check!" << std::endl;
0258         if (std::find(layerIdWithWireHistos.begin(), layerIdWithWireHistos.end(), layerId) ==
0259             layerIdWithWireHistos.end())
0260           layerIdWithWireHistos.push_back(layerId);
0261       }
0262       // Check for noise
0263       int nspikes = 0;
0264       for (int ibin = 0; ibin < hist.GetNbinsX(); ibin++) {
0265         if (hist.GetBinContent(ibin + 1) > hist.GetMaximum() * 0.001)
0266           nspikes++;
0267       }
0268       if (nspikes > 50) {
0269         std::cout << "Found a lot of (>50) small spikes in layer " << layerId
0270                   << ". Please check if all wires are functioning as expected!" << std::endl;
0271         if (std::find(layerIdWithWireHistos.begin(), layerIdWithWireHistos.end(), layerId) ==
0272             layerIdWithWireHistos.end())
0273           layerIdWithWireHistos.push_back(layerId);
0274       }
0275       hLayerPeaks.Fill(theTPPeakMap[layerId]);
0276     }
0277   }
0278 }
0279 
0280 void DTT0Calibration::endJob() {
0281   std::cout << "Analyzed " << nevents << " events" << std::endl;
0282 
0283   DTT0* t0sWRTChamber = new DTT0();
0284 
0285   if (debug)
0286     cout << "[DTT0CalibrationPerLayer]Writing histos to file!" << endl;
0287 
0288   theFile.cd();
0289   //hT0SectorHisto->Write();
0290   hLayerPeaks.Write();
0291   for (const auto& wHisto : theHistoWireMap) {
0292     wHisto.second.Write();
0293   }
0294   for (const auto& lHisto : theHistoLayerMap) {
0295     lHisto.second.Write();
0296   }
0297 
0298   if (debug)
0299     cout << "[DTT0Calibration] Compute and store t0 and sigma per wire" << endl;
0300 
0301   // Calculate uncertainties per wire (counting experiment)
0302   for (auto& wiret0 : theAbsoluteT0PerWire) {
0303     auto& wireId = wiret0.first;
0304     if (nDigiPerWire[wireId] > 1)
0305       theSigmaT0PerWire[wireId] = qK[wireId] / (nDigiPerWire[wireId] - 1);
0306     else
0307       theSigmaT0PerWire[wireId] = 999.;  // Only one measurement: uncertainty -> infinity
0308     // syst uncert
0309     //theSigmaT0PerWire[wireId] += pow(0.5, 2));
0310     // Every time the same measurement. Use Laplace estimator as estimation how propable it is to measure another value due to limited size of sample
0311     if (theSigmaT0PerWire[wireId] == 0) {
0312       theSigmaT0PerWire[wireId] += pow(1. / (nDigiPerWire[wireId] + 1), 2);
0313     }
0314   }
0315 
0316   // function to calculate unweighted means
0317   auto unweighted_mean_function = [](const std::list<double>& values, const std::list<double>& sigmas) {
0318     double mean = 0;
0319     for (auto& value : values) {
0320       mean += value;
0321     }
0322     mean /= values.size();
0323 
0324     double uncertainty = 0;
0325     for (auto& value : values) {
0326       uncertainty += pow(value - mean, 2);
0327     }
0328     uncertainty /= values.size();
0329     uncertainty = sqrt(uncertainty);
0330     return std::make_pair(mean, uncertainty);
0331   };
0332 
0333   // correct for odd-even effect in each super layer
0334   std::map<DTSuperLayerId, std::pair<double, double> > mean_sigma_even;
0335   std::map<DTSuperLayerId, std::pair<double, double> > mean_sigma_odd;
0336   for (const auto& superlayer : dtGeom->superLayers()) {
0337     const auto superlayer_id = superlayer->id();
0338     std::list<double> values_even;
0339     std::list<double> sigmas_even;
0340     std::list<double> values_odd;
0341     std::list<double> sigmas_odd;
0342 
0343     for (const auto& wiret0 : theAbsoluteT0PerWire) {
0344       const auto& wireId = wiret0.first;
0345       if (wireId.layerId().superlayerId() == superlayer_id) {
0346         const auto& t0 = wiret0.second / nDigiPerWire[wireId];
0347         if (wireId.layerId().layer() % 2) {
0348           values_odd.push_back(t0);
0349           sigmas_odd.push_back(sqrt(theSigmaT0PerWire[wireId]));
0350         } else {
0351           values_even.push_back(t0);
0352           sigmas_even.push_back(sqrt(theSigmaT0PerWire[wireId]));
0353         }
0354       }
0355     }
0356     // get mean and uncertainty
0357     mean_sigma_even.emplace(superlayer_id, unweighted_mean_function(values_even, sigmas_even));
0358     mean_sigma_odd.emplace(superlayer_id, unweighted_mean_function(values_odd, sigmas_odd));
0359   }
0360 
0361   // filter outliers
0362   for (const auto& superlayer : dtGeom->superLayers()) {
0363     const auto superlayer_id = superlayer->id();
0364     std::list<double> values_even;
0365     std::list<double> sigmas_even;
0366     std::list<double> values_odd;
0367     std::list<double> sigmas_odd;
0368 
0369     for (const auto& wiret0 : theAbsoluteT0PerWire) {
0370       const auto& wireId = wiret0.first;
0371       if (wireId.layerId().superlayerId() == superlayer_id) {
0372         const auto& t0 = wiret0.second / nDigiPerWire[wireId];
0373         if (wireId.layerId().layer() % 2 and
0374             abs(t0 - mean_sigma_odd[superlayer_id].first) < 2 * mean_sigma_odd[superlayer_id].second) {
0375           values_odd.push_back(t0);
0376           sigmas_odd.push_back(sqrt(theSigmaT0PerWire[wireId]));
0377         } else {
0378           if (abs(t0 - mean_sigma_even[superlayer_id].first) < 2 * mean_sigma_even[superlayer_id].second) {
0379             values_even.push_back(t0);
0380             sigmas_even.push_back(sqrt(theSigmaT0PerWire[wireId]));
0381           }
0382         }
0383       }
0384     }
0385     // get mean and uncertainty
0386     mean_sigma_even[superlayer_id] = unweighted_mean_function(values_even, sigmas_even);
0387     mean_sigma_odd[superlayer_id] = unweighted_mean_function(values_odd, sigmas_odd);
0388   }
0389 
0390   // apply correction
0391   for (auto& wiret0 : theAbsoluteT0PerWire) {
0392     const auto& wire_id = wiret0.first;
0393     const auto& superlayer_id = wiret0.first.layerId().superlayerId();
0394     const auto& layer = wiret0.first.layerId().layer();
0395     auto& t0 = wiret0.second;
0396     t0 /= nDigiPerWire[wire_id];
0397     if (not layer % 2)
0398       continue;
0399     // t0 is reference. Changing it changes the map
0400     t0 += mean_sigma_even[superlayer_id].first - mean_sigma_odd[superlayer_id].first;
0401     theSigmaT0PerWire[wire_id] +=
0402         pow(mean_sigma_odd[superlayer_id].second, 2) + pow(mean_sigma_even[superlayer_id].second, 2);
0403   }
0404 
0405   // get chamber mean
0406   std::map<DTChamberId, std::list<double> > values_per_chamber;
0407   std::map<DTChamberId, std::list<double> > sigmas_per_chamber;
0408   for (const auto& wire_t0 : theAbsoluteT0PerWire) {
0409     const auto& wire_id = wire_t0.first;
0410     const auto& chamber_id = wire_id.chamberId();
0411     const auto& t0 = wire_t0.second;
0412     values_per_chamber[chamber_id].push_back(t0);
0413     sigmas_per_chamber[chamber_id].push_back(sqrt(theSigmaT0PerWire[wire_id]));
0414   }
0415 
0416   std::map<DTChamberId, std::pair<double, double> > mean_per_chamber;
0417   for (const auto& chamber_mean : values_per_chamber) {
0418     const auto& chamber_id = chamber_mean.first;
0419     const auto& means = chamber_mean.second;
0420     const auto& sigmas = sigmas_per_chamber[chamber_id];
0421     mean_per_chamber.emplace(chamber_id, unweighted_mean_function(means, sigmas));
0422   }
0423 
0424   // calculate relative values
0425   for (const auto& wire_t0 : theAbsoluteT0PerWire) {
0426     const auto& wire_id = wire_t0.first;
0427     const auto& chamber_id = wire_id.chamberId();
0428     const auto& t0 = wire_t0.second;
0429     theRelativeT0PerWire.emplace(wire_id, t0 - mean_per_chamber[chamber_id].first);
0430     cout << "[DTT0Calibration] Wire " << wire_id << " has    t0 " << theRelativeT0PerWire[wire_id]
0431          << " (relative, after even-odd layer corrections)  "
0432          << "    sigma " << sqrt(theSigmaT0PerWire[wire_id]) << endl;
0433   }
0434 
0435   for (const auto& wire_t0 : theRelativeT0PerWire) {
0436     const auto& wire_id = wire_t0.first;
0437     const auto& t0 = wire_t0.second;
0438     t0sWRTChamber->set(wire_id, t0, sqrt(theSigmaT0PerWire[wire_id]), DTTimeUnits::counts);
0439   }
0440 
0441   ///Write the t0 map into DB
0442   if (debug)
0443     cout << "[DTT0Calibration]Writing values in DB!" << endl;
0444   // FIXME: to be read from cfg?
0445   string t0Record = "DTT0Rcd";
0446   // Write the t0 map to DB
0447   DTCalibDBUtils::writeToDB(t0Record, t0sWRTChamber);
0448   delete t0sWRTChamber;
0449 }
0450 
0451 string DTT0Calibration::getHistoName(const DTWireId& wId) const {
0452   string histoName;
0453   stringstream theStream;
0454   theStream << "Ch_" << wId.wheel() << "_" << wId.station() << "_" << wId.sector() << "_SL" << wId.superlayer() << "_L"
0455             << wId.layer() << "_W" << wId.wire() << "_hT0Histo";
0456   theStream >> histoName;
0457   return histoName;
0458 }
0459 
0460 string DTT0Calibration::getHistoName(const DTLayerId& lId) const {
0461   string histoName;
0462   stringstream theStream;
0463   theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector() << "_SL" << lId.superlayer() << "_L"
0464             << lId.layer() << "_hT0Histo";
0465   theStream >> histoName;
0466   return histoName;
0467 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!