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
SerialTaskQueue *	globalLuminosityBlocksQueue () final
SerialTaskQueue *	globalRunsQueue () final
bool	wantsGlobalLuminosityBlocks () const final
bool	wantsGlobalRuns () 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
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::vector< ModuleDescription const * > &modules, 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
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
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)
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<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)

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

{
   // 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 83 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  )

override

Fill the maps with t0 (by channel)

Perform the real analysis.

Definition at line 91 of file DTT0Calibration.cc.

References funct::abs(), DTSuperLayerId::chamberId(), gather_cfg::cout, debug, digiToken, dtGeom, eventsForLayerT0, eventsForWireT0, spr::find(), edm::EventSetup::get(), getHistoName(), create_public_lumi_plots::hist, hLayerPeaks, layerIdWithWireHistos, mK, mK_ref, nDigiPerWire, nDigiPerWire_ref, nevents, or, qK, rejectDigiFromPeak, DTChamberId::sector(), selSector, selWheel, jetUpdater_cfi::sort, spectrum, DTLayerId::superlayerId(), genVertex_cff::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)
      eventSetup.get<MuonGeometryRecord>().get(dtGeom);

   // 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 289 of file DTT0Calibration.cc.

References funct::abs(), DTTimeUnits::counts, gather_cfg::cout, debug, dtGeom, plotBeamSpotDB::first, hLayerPeaks, SiStripPI::mean, nDigiPerWire, nevents, funct::pow(), qK, edm::second(), DTT0::set(), mathSSE::sqrt(), DTGeometry::superLayers(), genVertex_cff::t0, theAbsoluteT0PerWire, theFile, theHistoLayerMap, theHistoWireMap, theRelativeT0PerWire, theSigmaT0PerWire, relativeConstraints::value, MuonErrorMatrixValues_cff::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 485 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 494 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 61 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 64 of file DTT0Calibration.h.

Referenced by analyze().

edm::ESHandle<DTGeometry> DTT0Calibration::dtGeom

private

Definition at line 126 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

unsigned int DTT0Calibration::eventsForLayerT0

private

Definition at line 72 of file DTT0Calibration.h.

Referenced by analyze().

unsigned int DTT0Calibration::eventsForWireT0

private

Definition at line 74 of file DTT0Calibration.h.

Referenced by analyze().

TH1D DTT0Calibration::hLayerPeaks

private

Definition at line 95 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

std::vector<DTLayerId> DTT0Calibration::layerIdWithWireHistos

private

Definition at line 101 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 109 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 110 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 107 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 108 of file DTT0Calibration.h.

Referenced by analyze().

unsigned int DTT0Calibration::nevents

private

Definition at line 70 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 111 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

unsigned int DTT0Calibration::rejectDigiFromPeak

private

Definition at line 83 of file DTT0Calibration.h.

Referenced by analyze().

int DTT0Calibration::selSector

private

Definition at line 89 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

int DTT0Calibration::selWheel

private

Definition at line 87 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

TSpectrum DTT0Calibration::spectrum

private

Definition at line 97 of file DTT0Calibration.h.

Referenced by analyze().

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

private

Definition at line 104 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

std::string DTT0Calibration::theCalibSector

private

Definition at line 88 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

std::string DTT0Calibration::theCalibWheel

private

Definition at line 86 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().

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

private

Definition at line 120 of file DTT0Calibration.h.

TFile DTT0Calibration::theFile

private

Definition at line 67 of file DTT0Calibration.h.

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

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

private

Definition at line 92 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 113 of file DTT0Calibration.h.

Referenced by analyze(), and endJob().

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

private

Definition at line 122 of file DTT0Calibration.h.

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

private

Definition at line 123 of file DTT0Calibration.h.

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

private

Definition at line 105 of file DTT0Calibration.h.

Referenced by endJob().

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

private

Definition at line 121 of file DTT0Calibration.h.

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

private

Definition at line 116 of file DTT0Calibration.h.

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

private

Definition at line 106 of file DTT0Calibration.h.

Referenced by endJob().

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

private

Definition at line 119 of file DTT0Calibration.h.

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

private

Definition at line 115 of file DTT0Calibration.h.

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

private

Definition at line 117 of file DTT0Calibration.h.

Referenced by analyze().

double DTT0Calibration::tpPeakWidth

private

Definition at line 77 of file DTT0Calibration.h.

Referenced by analyze().

double DTT0Calibration::tpPeakWidthPerLayer

private

Definition at line 80 of file DTT0Calibration.h.

Referenced by analyze().

std::vector<DTWireId> DTT0Calibration::wireIdWithHistos

private

Definition at line 100 of file DTT0Calibration.h.

Referenced by analyze(), and DTT0Calibration().