DTT0Calibration.cc

Go to the documentation of this file.

/*
 *  See header file for a description of this class.
 *
 *  $Date: 2012/05/11 17:17:17 $
 *  $Revision: 1.6 $
 *  \author S. Bolognesi - INFN Torino
 *  06/08/2008 Mofified by Antonio.Vilela.Pereira@cern.ch
 */

#include "CalibMuon/DTCalibration/plugins/DTT0Calibration.h"
#include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Geometry/Records/interface/MuonNumberingRecord.h"

#include "DataFormats/DTDigi/interface/DTDigiCollection.h"
#include "CondFormats/DTObjects/interface/DTT0.h"

#include <CondFormats/DTObjects/interface/DTTtrig.h>
#include <CondFormats/DataRecord/interface/DTTtrigRcd.h>

#include "TKey.h"
#include "TF1.h"

#include <cassert>

using namespace std;
using namespace edm;

// Constructor
DTT0Calibration::DTT0Calibration(const edm::ParameterSet& pset) :
   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));
   }
}

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

   theFile.Close();
}

void DTT0Calibration::analyze(const edm::Event & event, const edm::EventSetup& eventSetup) {
   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() {

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