#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 ( const edm::ParameterSet & pset )

Constructor.

Definition at line 34 of file DTT0Calibration.cc.

References gather_cfg::cout, debug, edm::ParameterSet::getUntrackedParameter(), phase1PixelTopology::layer, selSector, selWheel, relativeConstraints::station, theCalibSector, theCalibWheel, 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 ( )

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

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_, eventsForLayerT0, eventsForWireT0, spr::find(), edm::EventSetup::getHandle(), getHistoName(), gpuVertexFinder::hist, hLayerPeaks, layerIdWithWireHistos, mK, mK_ref, nDigiPerWire, nDigiPerWire_ref, nevents, or, qK, rejectDigiFromPeak, DTChamberId::sector(), selSector, selWheel, 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]);
     }
   }
 }

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(), FrontierCondition_GT_autoExpress_cfi::t0, theAbsoluteT0PerWire, theFile, theHistoLayerMap, theHistoWireMap, theRelativeT0PerWire, theSigmaT0PerWire, relativeConstraints::value, makeHLTPrescaleTable::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;
 }

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

private

Definition at line 451 of file DTT0Calibration.cc.

References 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;
 }

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

private

Definition at line 460 of file DTT0Calibration.cc.

References 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

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().

edm::EDGetTokenT<DTDigiCollection> DTT0Calibration::digiToken

private

Definition at line 62 of file DTT0Calibration.h.

Referenced by analyze().

edm::ESHandle<DTGeometry> DTT0Calibration::dtGeom

private

Definition at line 124 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 125 of file DTT0Calibration.h.

Referenced by analyze().

unsigned int DTT0Calibration::eventsForLayerT0

private

Definition at line 70 of file DTT0Calibration.h.

Referenced by analyze().

unsigned int DTT0Calibration::eventsForWireT0

private

Definition at line 72 of file DTT0Calibration.h.

Referenced by analyze().

TH1D DTT0Calibration::hLayerPeaks

private

Definition at line 93 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

std::vector<DTLayerId> DTT0Calibration::layerIdWithWireHistos

private

Definition at line 99 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 107 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 108 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 105 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 106 of file DTT0Calibration.h.

Referenced by analyze().

unsigned int DTT0Calibration::nevents

private

Definition at line 68 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 109 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

unsigned int DTT0Calibration::rejectDigiFromPeak

private

Definition at line 81 of file DTT0Calibration.h.

Referenced by analyze().

int DTT0Calibration::selSector

private

Definition at line 87 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

int DTT0Calibration::selWheel

private

Definition at line 85 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

TSpectrum DTT0Calibration::spectrum

private

Definition at line 95 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 102 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

std::string DTT0Calibration::theCalibSector

private

Definition at line 86 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

std::string DTT0Calibration::theCalibWheel

private

Definition at line 84 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

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

private

Definition at line 118 of file DTT0Calibration.h.

TFile DTT0Calibration::theFile

private

Definition at line 65 of file DTT0Calibration.h.

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

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

private

Definition at line 90 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 111 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 120 of file DTT0Calibration.h.

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

private

Definition at line 121 of file DTT0Calibration.h.

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

private

Definition at line 103 of file DTT0Calibration.h.

Referenced by endJob().

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

private

Definition at line 119 of file DTT0Calibration.h.

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

private

Definition at line 114 of file DTT0Calibration.h.

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

private

Definition at line 104 of file DTT0Calibration.h.

Referenced by endJob().

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

private

Definition at line 117 of file DTT0Calibration.h.

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

private

Definition at line 113 of file DTT0Calibration.h.

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

private

Definition at line 115 of file DTT0Calibration.h.

Referenced by analyze().

double DTT0Calibration::tpPeakWidth

private

Definition at line 75 of file DTT0Calibration.h.

Referenced by analyze().

double DTT0Calibration::tpPeakWidthPerLayer

private

Definition at line 78 of file DTT0Calibration.h.

Referenced by analyze().

std::vector<DTWireId> DTT0Calibration::wireIdWithHistos

private

Definition at line 98 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation