DTT0Calibration Class Reference

#include <DTT0Calibration.h>

Inheritance diagram for DTT0Calibration:

Public Member Functions
void	analyze (const edm::Event &event, const edm::EventSetup &eventSetup) override
	Fill the maps with t0 (by channel) More...
	DTT0Calibration (const edm::ParameterSet &pset)
	Constructor. More...
void	endJob () override
	Compute the mean and the RMS of the t0 from the maps and write them to the DB with channel granularity. More...
	~DTT0Calibration () override
	Destructor. More...
Public Member Functions inherited from edm::one::EDAnalyzer<>
	EDAnalyzer ()=default
	EDAnalyzer (const EDAnalyzer &)=delete
SerialTaskQueue *	globalLuminosityBlocksQueue () final
SerialTaskQueue *	globalRunsQueue () final
const EDAnalyzer &	operator= (const EDAnalyzer &)=delete
bool	wantsGlobalLuminosityBlocks () const final
bool	wantsGlobalRuns () const final
bool	wantsInputProcessBlocks () const final
bool	wantsProcessBlocks () const final
Public Member Functions inherited from edm::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
	EDAnalyzerBase ()
ModuleDescription const &	moduleDescription () const
bool	wantsStreamLuminosityBlocks () const
bool	wantsStreamRuns () const
	~EDAnalyzerBase () override
Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const
void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...
	EDConsumerBase ()
	EDConsumerBase (EDConsumerBase const &)=delete
	EDConsumerBase (EDConsumerBase &&)=default
ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const
std::vector< ESProxyIndex > const &	esGetTokenIndicesVector (edm::Transition iTrans) const
std::vector< ESRecordIndex > const &	esGetTokenRecordIndicesVector (edm::Transition iTrans) const
ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const
void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const
void	labelsForToken (EDGetToken iToken, Labels &oLabels) const
void	modulesWhoseProductsAreConsumed (std::array< std::vector< ModuleDescription const *> *, NumBranchTypes > &modulesAll, std::vector< ModuleProcessName > &modulesInPreviousProcesses, ProductRegistry const &preg, std::map< std::string, ModuleDescription const *> const &labelsToDesc, std::string const &processName) const
EDConsumerBase const &	operator= (EDConsumerBase const &)=delete
EDConsumerBase &	operator= (EDConsumerBase &&)=default
bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const
bool	registeredToConsumeMany (TypeID const &, BranchType) const
void	selectInputProcessBlocks (ProductRegistry const &productRegistry, ProcessBlockHelperBase const &processBlockHelperBase)
ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const
void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)
virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
std::string	getHistoName (const DTWireId &wId) const
std::string	getHistoName (const DTLayerId &lId) const

Private Attributes
bool	debug
edm::EDGetTokenT< DTDigiCollection >	digiToken
edm::ESHandle< DTGeometry >	dtGeom
const edm::ESGetToken< DTGeometry, MuonGeometryRecord >	dtGeomToken_
unsigned int	eventsForLayerT0
unsigned int	eventsForWireT0
TH1D	hLayerPeaks
std::vector< DTLayerId >	layerIdWithWireHistos
std::map< DTWireId, double >	mK
std::map< DTWireId, double >	mK_ref
std::map< DTWireId, int >	nDigiPerWire
std::map< DTWireId, int >	nDigiPerWire_ref
unsigned int	nevents
std::map< DTWireId, double >	qK
unsigned int	rejectDigiFromPeak
int	selSector
int	selWheel
TSpectrum	spectrum
std::map< DTWireId, double >	theAbsoluteT0PerWire
std::string	theCalibSector
std::string	theCalibWheel
std::map< DTChamberId, int >	theCountT0ByChamber
TFile	theFile
std::map< DTLayerId, TH1I >	theHistoLayerMap
std::map< DTWireId, TH1I >	theHistoWireMap
std::map< DTChamberId, double >	theMeanT0ByChamber
std::map< DTChamberId, double >	theRefT0ByChamber
std::map< DTWireId, double >	theRelativeT0PerWire
std::map< DTChamberId, double >	theSigmaT0ByChamber
std::map< std::string, double >	theSigmaT0LayerMap
std::map< DTWireId, double >	theSigmaT0PerWire
std::map< DTChamberId, double >	theSumT0ByChamber
std::map< std::string, double >	theT0LayerMap
std::map< DTLayerId, double >	theTPPeakMap
double	tpPeakWidth
double	tpPeakWidthPerLayer
std::vector< DTWireId >	wireIdWithHistos

Additional Inherited Members
Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase	ModuleType
Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels
Static Public Member Functions inherited from edm::one::EDAnalyzerBase
static const std::string &	baseType ()
static void	fillDescriptions (ConfigurationDescriptions &descriptions)
static void	prevalidate (ConfigurationDescriptions &descriptions)
Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)
template<BranchType B = InEvent>
EDConsumerBaseAdaptor< B >	consumes (edm::InputTag tag) noexcept
EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)
template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)
ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...
template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()
void	consumesMany (const TypeToGet &id)
template<BranchType B>
void	consumesMany (const TypeToGet &id)
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)
template<Transition Tr = Transition::Event>
constexpr auto	esConsumes ()
template<Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag tag)
template<Transition Tr = Transition::Event>
ESGetTokenGeneric	esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
	Used with EventSetupRecord::doGet. More...
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)
template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)
void	resetItemsToGetFrom (BranchType iType)

Detailed Description

Analyzer class computes the mean and RMS of t0 from pulses. Those values are written in the DB with cell granularity. The mean value for each channel is normalized to a reference time common to all the sector. The t0 of wires in odd layers are corrected for the relative difference between odd and even layers

Author

S. Bolognesi - INFN Torino

Definition at line 36 of file DTT0Calibration.h.

Constructor & Destructor Documentation

DTT0Calibration()

DTT0Calibration::DTT0Calibration

(

const edm::ParameterSet &

pset

)

Constructor.

Definition at line 34 of file DTT0Calibration.cc.

References gather_cfg::cout, debug, phase1PixelTopology::layer, muonDTDigis_cfi::pset, selSector, selWheel, relativeConstraints::station, theCalibSector, theCalibWheel, makeMuonMisalignmentScenario::wheel, and wireIdWithHistos.

    : debug(pset.getUntrackedParameter<bool>("debug")),
      digiToken(consumes<DTDigiCollection>(pset.getUntrackedParameter<string>("digiLabel"))),
      theFile(pset.getUntrackedParameter<string>("rootFileName", "DTT0PerLayer.root").c_str(), "RECREATE"),
      nevents(0),
      eventsForLayerT0(pset.getParameter<unsigned int>("eventsForLayerT0")),
      eventsForWireT0(pset.getParameter<unsigned int>("eventsForWireT0")),
      tpPeakWidth(pset.getParameter<double>("tpPeakWidth")),
      tpPeakWidthPerLayer(pset.getParameter<double>("tpPeakWidthPerLayer")),
      rejectDigiFromPeak(pset.getParameter<unsigned int>("rejectDigiFromPeak")),
      hLayerPeaks("hLayerPeaks", "", 3000, 0, 3000),
      spectrum(20),
      dtGeomToken_(esConsumes()) {
  // Get the debug parameter for verbose output
  if (debug)
    cout << "[DTT0Calibration]Constructor called!" << endl;

  theCalibWheel =
      pset.getUntrackedParameter<string>("calibWheel", "All");  //FIXME amke a vector of integer instead of a string
  if (theCalibWheel != "All") {
    stringstream linestr;
    int selWheel;
    linestr << theCalibWheel;
    linestr >> selWheel;
    cout << "[DTT0CalibrationPerLayer] chosen wheel " << selWheel << endl;
  }

  // Sector/s to calibrate
  theCalibSector =
      pset.getUntrackedParameter<string>("calibSector", "All");  //FIXME amke a vector of integer instead of a string
  if (theCalibSector != "All") {
    stringstream linestr;
    int selSector;
    linestr << theCalibSector;
    linestr >> selSector;
    cout << "[DTT0CalibrationPerLayer] chosen sector " << selSector << endl;
  }

  vector<string> defaultCell;
  const auto& cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell);
  for (const auto& cell : cellsWithHistos) {
    stringstream linestr;
    int wheel, sector, station, sl, layer, wire;
    linestr << cell;
    linestr >> wheel >> sector >> station >> sl >> layer >> wire;
    wireIdWithHistos.push_back(DTWireId(wheel, station, sector, sl, layer, wire));
  }
}

~DTT0Calibration()

DTT0Calibration::~DTT0Calibration ( )

override

Destructor.

Definition at line 84 of file DTT0Calibration.cc.

References gather_cfg::cout, debug, and theFile.

                                  {
  if (debug)
    cout << "[DTT0Calibration]Destructor called!" << endl;

  theFile.Close();
}

Member Function Documentation

analyze()

void DTT0Calibration::analyze ( const edm::Event &  event, const edm::EventSetup &  eventSetup  )

overridevirtual

Fill the maps with t0 (by channel)

Perform the real analysis.

Implements edm::one::EDAnalyzerBase.

Definition at line 92 of file DTT0Calibration.cc.

References funct::abs(), DTSuperLayerId::chamberId(), gather_cfg::cout, debug, digiToken, dtGeom, dtGeomToken_, options_cfi::eventSetup, eventsForLayerT0, eventsForWireT0, spr::find(), getHistoName(), gpuVertexFinder::hist, hLayerPeaks, layerIdWithWireHistos, mK, mK_ref, nDigiPerWire, nDigiPerWire_ref, nevents, or, qK, rejectDigiFromPeak, DTChamberId::sector(), selSector, selWheel, jetsAK4_CHS_cff::sort, spectrum, DTLayerId::superlayerId(), FrontierCondition_GT_autoExpress_cfi::t0, theAbsoluteT0PerWire, theCalibSector, theCalibWheel, theHistoLayerMap, theHistoWireMap, theTPPeakMap, tpPeakWidth, tpPeakWidthPerLayer, DTChamberId::wheel(), and wireIdWithHistos.

                                                                                    {
  nevents++;

  // Get the digis from the event
  Handle<DTDigiCollection> digis;
  event.getByToken(digiToken, digis);

  // Get the DT Geometry
  if (nevents == 1)
    dtGeom = eventSetup.getHandle(dtGeomToken_);

  // Iterate through all digi collections ordered by LayerId
  for (const auto& digis_per_layer : *digis) {
    //std::cout << __LINE__ << std::endl;
    // Get the iterators over the digis associated with this LayerId
    const DTDigiCollection::Range& digiRange = digis_per_layer.second;

    // Get the layerId
    const DTLayerId layerId = digis_per_layer.first;
    //const DTChamberId chamberId = layerId.superlayerId().chamberId();

    if ((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel))
      continue;
    if ((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector))
      continue;

    // Loop over all digis in the given layer
    for (DTDigiCollection::const_iterator digi = digiRange.first; digi != digiRange.second; ++digi) {
      const double t0 = digi->countsTDC();
      const DTWireId wireIdtmp(layerId, (*digi).wire());

      // Use first bunch of events to fill t0 per layer
      if (nevents <= eventsForLayerT0) {
        // If it doesn't exist, book it
        if (not theHistoLayerMap.count(layerId)) {
          theHistoLayerMap[layerId] = TH1I(getHistoName(layerId).c_str(),
                                           "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
                                           3000,
                                           0,
                                           3000);
          if (debug)
            cout << "  New T0 per Layer Histo: " << theHistoLayerMap[layerId].GetName() << endl;
        }
        theHistoLayerMap[layerId].Fill(t0);
      }

      // Use all the remaining events to compute t0 per wire
      if (nevents > eventsForLayerT0) {
        // Get the wireId
        const DTWireId wireId(layerId, (*digi).wire());
        if (debug) {
          cout << "   Wire: " << wireId << endl << "       time (TDC counts): " << (*digi).countsTDC() << endl;
        }

        //Fill the histos per wire for the chosen cells
        if (std::find(layerIdWithWireHistos.begin(), layerIdWithWireHistos.end(), layerId) !=
                layerIdWithWireHistos.end() or
            std::find(wireIdWithHistos.begin(), wireIdWithHistos.end(), wireId) != wireIdWithHistos.end()) {
          //If it doesn't exist, book it
          if (theHistoWireMap.count(wireId) == 0) {
            theHistoWireMap[wireId] = TH1I(getHistoName(wireId).c_str(),
                                           "T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",
                                           7000,
                                           0,
                                           7000);
            if (debug)
              cout << "  New T0 per wire Histo: " << theHistoWireMap[wireId].GetName() << endl;
          }
          theHistoWireMap[wireId].Fill(t0);
        }

        //Select per layer
        if (fabs(theTPPeakMap[layerId] - t0) > rejectDigiFromPeak) {
          if (debug)
            cout << "digi skipped because t0 too far from peak " << theTPPeakMap[layerId] << endl;
          continue;
        }

        //Use second bunch of events to compute a t0 reference per wire
        if (nevents <= (eventsForLayerT0 + eventsForWireT0)) {
          if (!nDigiPerWire_ref[wireId]) {
            mK_ref[wireId] = 0;
          }
          nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1;
          mK_ref[wireId] = mK_ref[wireId] + (t0 - mK_ref[wireId]) / nDigiPerWire_ref[wireId];
        }
        //Use last all the remaining events to compute the mean and sigma t0 per wire
        else if (nevents > (eventsForLayerT0 + eventsForWireT0)) {
          if (abs(t0 - mK_ref[wireId]) > tpPeakWidth)
            continue;
          if (!nDigiPerWire[wireId]) {
            theAbsoluteT0PerWire[wireId] = 0;
            qK[wireId] = 0;
            mK[wireId] = 0;
          }
          nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1;
          theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0;
          qK[wireId] =
              qK[wireId] + ((nDigiPerWire[wireId] - 1) * (t0 - mK[wireId]) * (t0 - mK[wireId]) / nDigiPerWire[wireId]);
          mK[wireId] = mK[wireId] + (t0 - mK[wireId]) / nDigiPerWire[wireId];
        }
      }  //end if(nevents>1000)
    }    //end loop on digi
  }      //end loop on layer

  //Use the t0 per layer histos to have an indication about the t0 position
  if (nevents == eventsForLayerT0) {
    for (const auto& lHisto : theHistoLayerMap) {
      const auto& layerId = lHisto.first;
      const auto& hist = lHisto.second;
      if (debug)
        cout << "Reading histogram " << hist.GetName() << " with mean " << hist.GetMean() << " and RMS "
             << hist.GetRMS() << endl;

      //Find peaks
      int npeaks = spectrum.Search(&hist, (tpPeakWidthPerLayer / 2.), "", 0.3);

      double* peaks = spectrum.GetPositionX();
      //Put in a std::vector<float>
      vector<double> peakMeans(peaks, peaks + npeaks);
      //Sort the peaks in ascending order
      sort(peakMeans.begin(), peakMeans.end());

      if (peakMeans.empty()) {
        theTPPeakMap[layerId] = hist.GetMaximumBin();
        std::cout << "No peaks found by peakfinder in layer " << layerId << ". Taking maximum bin at "
                  << theTPPeakMap[layerId] << ". Please check!" << std::endl;
        layerIdWithWireHistos.push_back(layerId);
      } else if (fabs(hist.GetXaxis()->FindBin(peakMeans.front()) - hist.GetXaxis()->FindBin(peakMeans.back())) <
                 rejectDigiFromPeak) {
        theTPPeakMap[layerId] = peakMeans[peakMeans.size() / 2];
      } else {
        bool peak_set = false;
        for (const auto& peak : peakMeans) {
          // Skip if at low edge
          if (peak - tpPeakWidthPerLayer <= 0)
            continue;
          // Get integral of peak
          double sum = 0;
          for (int ibin = peak - tpPeakWidthPerLayer; ibin < peak + tpPeakWidthPerLayer; ibin++) {
            sum += hist.GetBinContent(ibin);
          }
          // Skip if peak too small
          if (sum < hist.GetMaximum() / 2)
            continue;

          // Passed all cuts
          theTPPeakMap[layerId] = peak;
          peak_set = true;
          break;
        }
        if (peak_set) {
          std::cout << "Peaks to far away from each other in layer " << layerId
                    << ". Maybe cross talk? Taking first good peak at " << theTPPeakMap[layerId] << ". Please check!"
                    << std::endl;
          layerIdWithWireHistos.push_back(layerId);
        } else {
          theTPPeakMap[layerId] = hist.GetMaximumBin();
          std::cout << "Peaks to far away from each other in layer " << layerId
                    << " and no good peak found. Taking maximum bin at " << theTPPeakMap[layerId] << ". Please check!"
                    << std::endl;
          layerIdWithWireHistos.push_back(layerId);
        }
      }
      if (peakMeans.size() > 5) {
        std::cout << "Found more than 5 peaks in layer " << layerId << ". Please check!" << std::endl;
        if (std::find(layerIdWithWireHistos.begin(), layerIdWithWireHistos.end(), layerId) ==
            layerIdWithWireHistos.end())
          layerIdWithWireHistos.push_back(layerId);
      }
      // Check for noise
      int nspikes = 0;
      for (int ibin = 0; ibin < hist.GetNbinsX(); ibin++) {
        if (hist.GetBinContent(ibin + 1) > hist.GetMaximum() * 0.001)
          nspikes++;
      }
      if (nspikes > 50) {
        std::cout << "Found a lot of (>50) small spikes in layer " << layerId
                  << ". Please check if all wires are functioning as expected!" << std::endl;
        if (std::find(layerIdWithWireHistos.begin(), layerIdWithWireHistos.end(), layerId) ==
            layerIdWithWireHistos.end())
          layerIdWithWireHistos.push_back(layerId);
      }
      hLayerPeaks.Fill(theTPPeakMap[layerId]);
    }
  }
}

endJob()

void DTT0Calibration::endJob ( void  )

overridevirtual

Compute the mean and the RMS of the t0 from the maps and write them to the DB with channel granularity.

Write the t0 map into DB

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 280 of file DTT0Calibration.cc.

References funct::abs(), DTTimeUnits::counts, gather_cfg::cout, debug, dtGeom, first, hLayerPeaks, phase1PixelTopology::layer, SiStripPI::mean, compareTotals::means, nDigiPerWire, nevents, funct::pow(), qK, edm::second(), DTT0::set(), mathSSE::sqrt(), DTGeometry::superLayers(), FrontierCondition_GT_autoExpress_cfi::t0, theAbsoluteT0PerWire, theFile, theHistoLayerMap, theHistoWireMap, theRelativeT0PerWire, theSigmaT0PerWire, relativeConstraints::value, contentValuesCheck::values, and DTCalibDBUtils::writeToDB().

Referenced by o2olib.O2ORunMgr::executeJob().

                             {
  std::cout << "Analyzed " << nevents << " events" << std::endl;

  DTT0* t0sWRTChamber = new DTT0();

  if (debug)
    cout << "[DTT0CalibrationPerLayer]Writing histos to file!" << endl;

  theFile.cd();
  //hT0SectorHisto->Write();
  hLayerPeaks.Write();
  for (const auto& wHisto : theHistoWireMap) {
    wHisto.second.Write();
  }
  for (const auto& lHisto : theHistoLayerMap) {
    lHisto.second.Write();
  }

  if (debug)
    cout << "[DTT0Calibration] Compute and store t0 and sigma per wire" << endl;

  // Calculate uncertainties per wire (counting experiment)
  for (auto& wiret0 : theAbsoluteT0PerWire) {
    auto& wireId = wiret0.first;
    if (nDigiPerWire[wireId] > 1)
      theSigmaT0PerWire[wireId] = qK[wireId] / (nDigiPerWire[wireId] - 1);
    else
      theSigmaT0PerWire[wireId] = 999.;  // Only one measurement: uncertainty -> infinity
    // syst uncert
    //theSigmaT0PerWire[wireId] += pow(0.5, 2));
    // Every time the same measurement. Use Laplace estimator as estimation how propable it is to measure another value due to limited size of sample
    if (theSigmaT0PerWire[wireId] == 0) {
      theSigmaT0PerWire[wireId] += pow(1. / (nDigiPerWire[wireId] + 1), 2);
    }
  }

  // function to calculate unweighted means
  auto unweighted_mean_function = [](const std::list<double>& values, const std::list<double>& sigmas) {
    double mean = 0;
    for (auto& value : values) {
      mean += value;
    }
    mean /= values.size();

    double uncertainty = 0;
    for (auto& value : values) {
      uncertainty += pow(value - mean, 2);
    }
    uncertainty /= values.size();
    uncertainty = sqrt(uncertainty);
    return std::make_pair(mean, uncertainty);
  };

  // correct for odd-even effect in each super layer
  std::map<DTSuperLayerId, std::pair<double, double> > mean_sigma_even;
  std::map<DTSuperLayerId, std::pair<double, double> > mean_sigma_odd;
  for (const auto& superlayer : dtGeom->superLayers()) {
    const auto superlayer_id = superlayer->id();
    std::list<double> values_even;
    std::list<double> sigmas_even;
    std::list<double> values_odd;
    std::list<double> sigmas_odd;

    for (const auto& wiret0 : theAbsoluteT0PerWire) {
      const auto& wireId = wiret0.first;
      if (wireId.layerId().superlayerId() == superlayer_id) {
        const auto& t0 = wiret0.second / nDigiPerWire[wireId];
        if (wireId.layerId().layer() % 2) {
          values_odd.push_back(t0);
          sigmas_odd.push_back(sqrt(theSigmaT0PerWire[wireId]));
        } else {
          values_even.push_back(t0);
          sigmas_even.push_back(sqrt(theSigmaT0PerWire[wireId]));
        }
      }
    }
    // get mean and uncertainty
    mean_sigma_even.emplace(superlayer_id, unweighted_mean_function(values_even, sigmas_even));
    mean_sigma_odd.emplace(superlayer_id, unweighted_mean_function(values_odd, sigmas_odd));
  }

  // filter outliers
  for (const auto& superlayer : dtGeom->superLayers()) {
    const auto superlayer_id = superlayer->id();
    std::list<double> values_even;
    std::list<double> sigmas_even;
    std::list<double> values_odd;
    std::list<double> sigmas_odd;

    for (const auto& wiret0 : theAbsoluteT0PerWire) {
      const auto& wireId = wiret0.first;
      if (wireId.layerId().superlayerId() == superlayer_id) {
        const auto& t0 = wiret0.second / nDigiPerWire[wireId];
        if (wireId.layerId().layer() % 2 and
            abs(t0 - mean_sigma_odd[superlayer_id].first) < 2 * mean_sigma_odd[superlayer_id].second) {
          values_odd.push_back(t0);
          sigmas_odd.push_back(sqrt(theSigmaT0PerWire[wireId]));
        } else {
          if (abs(t0 - mean_sigma_even[superlayer_id].first) < 2 * mean_sigma_even[superlayer_id].second) {
            values_even.push_back(t0);
            sigmas_even.push_back(sqrt(theSigmaT0PerWire[wireId]));
          }
        }
      }
    }
    // get mean and uncertainty
    mean_sigma_even[superlayer_id] = unweighted_mean_function(values_even, sigmas_even);
    mean_sigma_odd[superlayer_id] = unweighted_mean_function(values_odd, sigmas_odd);
  }

  // apply correction
  for (auto& wiret0 : theAbsoluteT0PerWire) {
    const auto& wire_id = wiret0.first;
    const auto& superlayer_id = wiret0.first.layerId().superlayerId();
    const auto& layer = wiret0.first.layerId().layer();
    auto& t0 = wiret0.second;
    t0 /= nDigiPerWire[wire_id];
    if (not layer % 2)
      continue;
    // t0 is reference. Changing it changes the map
    t0 += mean_sigma_even[superlayer_id].first - mean_sigma_odd[superlayer_id].first;
    theSigmaT0PerWire[wire_id] +=
        pow(mean_sigma_odd[superlayer_id].second, 2) + pow(mean_sigma_even[superlayer_id].second, 2);
  }

  // get chamber mean
  std::map<DTChamberId, std::list<double> > values_per_chamber;
  std::map<DTChamberId, std::list<double> > sigmas_per_chamber;
  for (const auto& wire_t0 : theAbsoluteT0PerWire) {
    const auto& wire_id = wire_t0.first;
    const auto& chamber_id = wire_id.chamberId();
    const auto& t0 = wire_t0.second;
    values_per_chamber[chamber_id].push_back(t0);
    sigmas_per_chamber[chamber_id].push_back(sqrt(theSigmaT0PerWire[wire_id]));
  }

  std::map<DTChamberId, std::pair<double, double> > mean_per_chamber;
  for (const auto& chamber_mean : values_per_chamber) {
    const auto& chamber_id = chamber_mean.first;
    const auto& means = chamber_mean.second;
    const auto& sigmas = sigmas_per_chamber[chamber_id];
    mean_per_chamber.emplace(chamber_id, unweighted_mean_function(means, sigmas));
  }

  // calculate relative values
  for (const auto& wire_t0 : theAbsoluteT0PerWire) {
    const auto& wire_id = wire_t0.first;
    const auto& chamber_id = wire_id.chamberId();
    const auto& t0 = wire_t0.second;
    theRelativeT0PerWire.emplace(wire_id, t0 - mean_per_chamber[chamber_id].first);
    cout << "[DTT0Calibration] Wire " << wire_id << " has    t0 " << theRelativeT0PerWire[wire_id]
         << " (relative, after even-odd layer corrections)  "
         << "    sigma " << sqrt(theSigmaT0PerWire[wire_id]) << endl;
  }

  for (const auto& wire_t0 : theRelativeT0PerWire) {
    const auto& wire_id = wire_t0.first;
    const auto& t0 = wire_t0.second;
    t0sWRTChamber->set(wire_id, t0, sqrt(theSigmaT0PerWire[wire_id]), DTTimeUnits::counts);
  }

  if (debug)
    cout << "[DTT0Calibration]Writing values in DB!" << endl;
  // FIXME: to be read from cfg?
  string t0Record = "DTT0Rcd";
  // Write the t0 map to DB
  DTCalibDBUtils::writeToDB(t0Record, t0sWRTChamber);
  delete t0sWRTChamber;
}

getHistoName() [1/2]

string DTT0Calibration::getHistoName ( const DTWireId &  wId ) const

private

Definition at line 451 of file DTT0Calibration.cc.

References HltBtagPostValidation_cff::histoName, DTLayerId::layer(), DTChamberId::sector(), DTChamberId::station(), DTSuperLayerId::superlayer(), DTChamberId::wheel(), and DTWireId::wire().

Referenced by analyze().

                                                              {
  string histoName;
  stringstream theStream;
  theStream << "Ch_" << wId.wheel() << "_" << wId.station() << "_" << wId.sector() << "_SL" << wId.superlayer() << "_L"
            << wId.layer() << "_W" << wId.wire() << "_hT0Histo";
  theStream >> histoName;
  return histoName;
}

getHistoName() [2/2]

string DTT0Calibration::getHistoName ( const DTLayerId &  lId ) const

private

Definition at line 460 of file DTT0Calibration.cc.

References HltBtagPostValidation_cff::histoName, DTLayerId::layer(), DTChamberId::sector(), DTChamberId::station(), DTSuperLayerId::superlayer(), and DTChamberId::wheel().

                                                               {
  string histoName;
  stringstream theStream;
  theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector() << "_SL" << lId.superlayer() << "_L"
            << lId.layer() << "_hT0Histo";
  theStream >> histoName;
  return histoName;
}

Member Data Documentation

debug

bool DTT0Calibration::debug

private

Definition at line 59 of file DTT0Calibration.h.

Referenced by analyze(), rrapi.RRApi::dprint(), util.rrapi.RRApi::dprint(), DTT0Calibration(), endJob(), pkg.AbstractPkg::generate(), rrapi.RRApi::get(), util.rrapi.RRApi::get(), pkg.AbstractPkg::get_kwds(), runTauIdMVA.TauIDEmbedder::loadMVA_WPs_run2_2017(), runTauIdMVA.TauIDEmbedder::runTauID(), crabFunctions.CrabController::submit(), pkg.AbstractPkg::write(), and ~DTT0Calibration().

digiToken

edm::EDGetTokenT<DTDigiCollection> DTT0Calibration::digiToken

private

Definition at line 62 of file DTT0Calibration.h.

Referenced by analyze().

dtGeom

edm::ESHandle<DTGeometry> DTT0Calibration::dtGeom

private

Definition at line 124 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

dtGeomToken_

const edm::ESGetToken<DTGeometry, MuonGeometryRecord> DTT0Calibration::dtGeomToken_

private

Definition at line 125 of file DTT0Calibration.h.

Referenced by analyze().

eventsForLayerT0

unsigned int DTT0Calibration::eventsForLayerT0

private

Definition at line 70 of file DTT0Calibration.h.

Referenced by analyze().

eventsForWireT0

unsigned int DTT0Calibration::eventsForWireT0

private

Definition at line 72 of file DTT0Calibration.h.

Referenced by analyze().

hLayerPeaks

TH1D DTT0Calibration::hLayerPeaks

private

Definition at line 93 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

layerIdWithWireHistos

std::vector<DTLayerId> DTT0Calibration::layerIdWithWireHistos

private

Definition at line 99 of file DTT0Calibration.h.

Referenced by analyze().

mK

std::map<DTWireId, double> DTT0Calibration::mK

private

Definition at line 107 of file DTT0Calibration.h.

Referenced by analyze().

mK_ref

std::map<DTWireId, double> DTT0Calibration::mK_ref

private

Definition at line 108 of file DTT0Calibration.h.

Referenced by analyze().

nDigiPerWire

std::map<DTWireId, int> DTT0Calibration::nDigiPerWire

private

Definition at line 105 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

nDigiPerWire_ref

std::map<DTWireId, int> DTT0Calibration::nDigiPerWire_ref

private

Definition at line 106 of file DTT0Calibration.h.

Referenced by analyze().

nevents

unsigned int DTT0Calibration::nevents

private

Definition at line 68 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

qK

std::map<DTWireId, double> DTT0Calibration::qK

private

Definition at line 109 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

rejectDigiFromPeak

unsigned int DTT0Calibration::rejectDigiFromPeak

private

Definition at line 81 of file DTT0Calibration.h.

Referenced by analyze().

selSector

int DTT0Calibration::selSector

private

Definition at line 87 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

selWheel

int DTT0Calibration::selWheel

private

Definition at line 85 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

spectrum

TSpectrum DTT0Calibration::spectrum

private

Definition at line 95 of file DTT0Calibration.h.

Referenced by analyze().

theAbsoluteT0PerWire

std::map<DTWireId, double> DTT0Calibration::theAbsoluteT0PerWire

private

Definition at line 102 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

theCalibSector

std::string DTT0Calibration::theCalibSector

private

Definition at line 86 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

theCalibWheel

std::string DTT0Calibration::theCalibWheel

private

Definition at line 84 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

theCountT0ByChamber

std::map<DTChamberId, int> DTT0Calibration::theCountT0ByChamber

private

Definition at line 118 of file DTT0Calibration.h.

theFile

TFile DTT0Calibration::theFile

private

Definition at line 65 of file DTT0Calibration.h.

Referenced by endJob(), and ~DTT0Calibration().

theHistoLayerMap

std::map<DTLayerId, TH1I> DTT0Calibration::theHistoLayerMap

private

Definition at line 90 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

theHistoWireMap

std::map<DTWireId, TH1I> DTT0Calibration::theHistoWireMap

private

Definition at line 111 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

theMeanT0ByChamber

std::map<DTChamberId, double> DTT0Calibration::theMeanT0ByChamber

private

Definition at line 120 of file DTT0Calibration.h.

theRefT0ByChamber

std::map<DTChamberId, double> DTT0Calibration::theRefT0ByChamber

private

Definition at line 121 of file DTT0Calibration.h.

theRelativeT0PerWire

std::map<DTWireId, double> DTT0Calibration::theRelativeT0PerWire

private

Definition at line 103 of file DTT0Calibration.h.

Referenced by endJob().

theSigmaT0ByChamber

std::map<DTChamberId, double> DTT0Calibration::theSigmaT0ByChamber

private

Definition at line 119 of file DTT0Calibration.h.

theSigmaT0LayerMap

std::map<std::string, double> DTT0Calibration::theSigmaT0LayerMap

private

Definition at line 114 of file DTT0Calibration.h.

theSigmaT0PerWire

std::map<DTWireId, double> DTT0Calibration::theSigmaT0PerWire

private

Definition at line 104 of file DTT0Calibration.h.

Referenced by endJob().

theSumT0ByChamber

std::map<DTChamberId, double> DTT0Calibration::theSumT0ByChamber

private

Definition at line 117 of file DTT0Calibration.h.

theT0LayerMap

std::map<std::string, double> DTT0Calibration::theT0LayerMap

private

Definition at line 113 of file DTT0Calibration.h.

theTPPeakMap

std::map<DTLayerId, double> DTT0Calibration::theTPPeakMap

private

Definition at line 115 of file DTT0Calibration.h.

Referenced by analyze().

tpPeakWidth

double DTT0Calibration::tpPeakWidth

private

Definition at line 75 of file DTT0Calibration.h.

Referenced by analyze().

tpPeakWidthPerLayer

double DTT0Calibration::tpPeakWidthPerLayer

private

Definition at line 78 of file DTT0Calibration.h.

Referenced by analyze().

wireIdWithHistos

std::vector<DTWireId> DTT0Calibration::wireIdWithHistos

private

Definition at line 98 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().