DTT0CalibrationRMS.cc

Go to the documentation of this file.

/*
 *  See header file for a description of this class.
 *
 *  \author S. Bolognesi - INFN Torino
 */
#include "CalibMuon/DTCalibration/plugins/DTT0CalibrationRMS.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 "TH1I.h"
#include "TFile.h"
#include "TKey.h"

using namespace std;
using namespace edm;
// using namespace cond;

// Constructor
DTT0CalibrationRMS::DTT0CalibrationRMS(const edm::ParameterSet& pset) {
  // Get the debug parameter for verbose output
  debug = pset.getUntrackedParameter<bool>("debug");
  if(debug)
    cout << "[DTT0CalibrationRMS]Constructor called!" << endl;

  // Get the label to retrieve digis from the event
  digiLabel = pset.getUntrackedParameter<string>("digiLabel");

  // The root file which contain the histos per layer
  string rootFileName = pset.getUntrackedParameter<string>("rootFileName","DTT0PerLayer.root");
  theFile = new TFile(rootFileName.c_str(), "RECREATE");

  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 << "[DTT0CalibrationRMSPerLayer] 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 << "[DTT0CalibrationRMSPerLayer] chosen sector " << selSector << endl;
  }

  vector<string> defaultCell;
  defaultCell.push_back("None");
  cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell);
  for(vector<string>::const_iterator cell = cellsWithHistos.begin(); cell != cellsWithHistos.end(); ++cell){
    if((*cell) != "None"){
      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));
    }
  }

  hT0SectorHisto=nullptr;

  nevents=0;
  eventsForLayerT0 = pset.getParameter<unsigned int>("eventsForLayerT0");
  eventsForWireT0 = pset.getParameter<unsigned int>("eventsForWireT0");
  rejectDigiFromPeak = pset.getParameter<unsigned int>("rejectDigiFromPeak");
  tpPeakWidth = pset.getParameter<double>("tpPeakWidth");
  //useReferenceWireInLayer_ = true;
  correctByChamberMean_ = pset.getParameter<bool>("correctByChamberMean");
}

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

  theFile->Close();
}

void DTT0CalibrationRMS::analyze(const edm::Event & event, const edm::EventSetup& eventSetup) {
  if(debug || event.id().event() % 500==0)
    cout << "--- [DTT0CalibrationRMS] Analysing Event: #Run: " << event.id().run()
     << " #Event: " << event.id().event() << endl;
  nevents++;

  // Get the digis from the event
  Handle<DTDigiCollection> digis; 
  event.getByLabel(digiLabel, digis);

  // Get the DT Geometry
  eventSetup.get<MuonGeometryRecord>().get(dtGeom);

  // Iterate through all digi collections ordered by LayerId   
  DTDigiCollection::DigiRangeIterator dtLayerIt;
  for (dtLayerIt = digis->begin();
       dtLayerIt != digis->end();
       ++dtLayerIt){
    // Get the iterators over the digis associated with this LayerId
    const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
  
    // Get the layerId
    const DTLayerId layerId = (*dtLayerIt).first; //FIXME: check to be in the right sector

    if((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel))
      continue;
    if((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector))
      continue;
 
    //if(debug) {
    //  cout << "Layer " << layerId<<" with "<<distance(digiRange.first, digiRange.second)<<" digi"<<endl;
    //}

    // Loop over all digis in the given layer
    for (DTDigiCollection::const_iterator digi = digiRange.first;
     digi != digiRange.second;
     ++digi) {
      double t0 = (*digi).countsTDC();

      //Use first bunch of events to fill t0 per layer
      if(nevents < eventsForLayerT0){
    //Get the per-layer histo from the map
    TH1I *hT0LayerHisto = theHistoLayerMap[layerId];
    //If it doesn't exist, book it
    if(hT0LayerHisto == nullptr){
      theFile->cd();
      hT0LayerHisto = new TH1I(getHistoName(layerId).c_str(),
                   "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
                   200, t0-100, t0+100);
      if(debug)
        cout << "  New T0 per Layer Histo: " << hT0LayerHisto->GetName() << endl;
      theHistoLayerMap[layerId] = hT0LayerHisto;
    }
    
    //Fill the histos
    theFile->cd();
    if(hT0LayerHisto != nullptr) {
      //  if(debug)
      // cout<<"Filling histo "<<hT0LayerHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl;
      hT0LayerHisto->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
    vector<DTWireId>::iterator it_wire = find(wireIdWithHistos.begin(),wireIdWithHistos.end(),wireId);
    if(it_wire != wireIdWithHistos.end()){
          if(theHistoWireMap.find(wireId) == theHistoWireMap.end()){
            theHistoWireMap[wireId] = new 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;
          }
          if(theHistoWireMap_ref.find(wireId) == theHistoWireMap_ref.end()){
            theHistoWireMap_ref[wireId] = new TH1I((getHistoName(wireId) + "_ref").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_ref[wireId])->GetName() << endl;
          }

          TH1I* hT0WireHisto = theHistoWireMap[wireId];
      //Fill the histos
      theFile->cd();
      if(hT0WireHisto) hT0WireHisto->Fill(t0);
    }

    //Check the tzero has reasonable value
    if(abs(hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin()) - t0) > rejectDigiFromPeak){
      if(debug)
        cout<<"digi skipped because t0 per sector "<<hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;
      continue;
    }

    //Use second bunch of events to compute a t0 reference per wire
    if(nevents< (eventsForLayerT0 + eventsForWireT0)){
          //Fill reference wire histos
          if(it_wire != wireIdWithHistos.end()){
            TH1I* hT0WireHisto_ref = theHistoWireMap_ref[wireId];
            theFile->cd();
            if(hT0WireHisto_ref) hT0WireHisto_ref->Fill(t0); 
          } 
      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;
      //theSigmaT0PerWire[wireId] = theSigmaT0PerWire[wireId] + (t0*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(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
    lHisto != theHistoLayerMap.end();
    ++lHisto){
      if(debug)
    cout<<"Reading histogram "<<(*lHisto).second->GetName()<<" with mean "<<(*lHisto).second->GetMean()<<" and RMS "<<(*lHisto).second->GetRMS();

      //Take the mean as a first t0 estimation
      if((*lHisto).second->GetRMS()<5.0){
    if(hT0SectorHisto == nullptr){
      hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 
                    //20, (*lHisto).second->GetMean()-100, (*lHisto).second->GetMean()+100);
                    700, 0, 7000);
    }
    if(debug)
      cout<<" accepted"<<endl;
    hT0SectorHisto->Fill((*lHisto).second->GetMean());
      }
      //Take the mean of noise + 400ns as a first t0 estimation
      // if((*lHisto).second->GetRMS()>10.0 && ((*lHisto).second->GetRMS()<15.0)){
//  double t0_estim = (*lHisto).second->GetMean() + 400;
//  if(hT0SectorHisto == 0){
//    hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 
//                  //20, t0_estim-100, t0_estim+100);
//                  700, 0, 7000);
//  }
//  if(debug)
//    cout<<" accepted + 400ns"<<endl;
//  hT0SectorHisto->Fill((*lHisto).second->GetMean() + 400);
//       }
      if(debug)
    cout<<endl;

      theT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetMean();
      theSigmaT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetRMS();
    }
    if(!hT0SectorHisto){
      cout<<"[DTT0CalibrationRMS]: All the t0 per layer are still uncorrect: trying with greater number of events"<<endl;
      eventsForLayerT0 = eventsForLayerT0*2;
      return;
    }
    if(debug)
      cout<<"[DTT0CalibrationRMS] t0 reference for this sector "<<
    hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;
  } 
}


void DTT0CalibrationRMS::endJob() {

  DTT0* t0sAbsolute = new DTT0();
  DTT0* t0sRelative = new DTT0();
  DTT0* t0sWRTChamber = new DTT0();

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

  theFile->cd();
  theFile->WriteTObject(hT0SectorHisto);
  //hT0SectorHisto->Write();
  for(map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap.begin();
      wHisto != theHistoWireMap.end();
      ++wHisto) {
    (*wHisto).second->Write(); 
  }
  for(map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap_ref.begin();
      wHisto != theHistoWireMap_ref.end();
      ++wHisto) {
    (*wHisto).second->Write();
  }
  for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
      lHisto != theHistoLayerMap.end();
      ++lHisto) {
    (*lHisto).second->Write(); 
  }  

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

  for(map<DTWireId, double>::const_iterator wiret0 = theAbsoluteT0PerWire.begin();
      wiret0 != theAbsoluteT0PerWire.end();
      ++wiret0){
    if(nDigiPerWire[(*wiret0).first]){
      double t0 = (*wiret0).second/nDigiPerWire[(*wiret0).first];

      theRelativeT0PerWire[(*wiret0).first] = t0 - hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin());

      //theSigmaT0PerWire[(*wiret0).first] = sqrt((theSigmaT0PerWire[(*wiret0).first] / nDigiPerWire[(*wiret0).first]) - t0*t0);
      theSigmaT0PerWire[(*wiret0).first] = sqrt(qK[(*wiret0).first]/nDigiPerWire[(*wiret0).first]);

      cout << "Wire " << (*wiret0).first << " has t0 " << t0 << "(absolute) "
                                                       << theRelativeT0PerWire[(*wiret0).first] << "(relative)"
                                         << "    sigma " << theSigmaT0PerWire[(*wiret0).first] << endl;

      t0sAbsolute->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first],DTTimeUnits::counts); 
    }
    else{
      cout<<"[DTT0CalibrationRMS] ERROR: no digis in wire "<<(*wiret0).first<<endl;
      abort();
    }
  }

  if(correctByChamberMean_){
     // Get all the sls from the setup
     const vector<const DTSuperLayer*> superLayers = dtGeom->superLayers();     
     // Loop over all SLs
     for(vector<const DTSuperLayer*>::const_iterator  sl = superLayers.begin();
       sl != superLayers.end(); sl++) {


    //Compute mean for odd and even superlayers
    double oddLayersMean=0;
    double evenLayersMean=0; 
    double oddLayersDen=0;
    double evenLayersDen=0;
    for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
          wiret0 != theRelativeT0PerWire.end();
          ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
          if(debug)
         cout<<"[DTT0CalibrationRMS] Superlayer "<<(*sl)->id()
            <<"layer " <<(*wiret0).first.layerId().layer()<<" with "<<(*wiret0).second<<endl;
          if(((*wiret0).first.layerId().layer()) % 2){
         oddLayersMean = oddLayersMean + (*wiret0).second;
         oddLayersDen++;
          }
          else{
         evenLayersMean = evenLayersMean + (*wiret0).second;
         evenLayersDen++;
          }
       }
    }
    oddLayersMean = oddLayersMean/oddLayersDen;
    evenLayersMean = evenLayersMean/evenLayersDen;
    //if(debug && oddLayersMean)
    cout<<"[DTT0CalibrationRMS] Relative T0 mean for  odd layers "<<oddLayersMean<<"  even layers"<<evenLayersMean<<endl;

    //Compute sigma for odd and even superlayers
    double oddLayersSigma=0;
    double evenLayersSigma=0;
    for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
          wiret0 != theRelativeT0PerWire.end();
          ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
          if(((*wiret0).first.layerId().layer()) % 2){
         oddLayersSigma = oddLayersSigma + ((*wiret0).second - oddLayersMean) * ((*wiret0).second - oddLayersMean);
          }
          else{
         evenLayersSigma = evenLayersSigma + ((*wiret0).second - evenLayersMean) * ((*wiret0).second - evenLayersMean);
          }
       }
    }
    oddLayersSigma = oddLayersSigma/oddLayersDen;
    evenLayersSigma = evenLayersSigma/evenLayersDen;
    oddLayersSigma = sqrt(oddLayersSigma);
    evenLayersSigma = sqrt(evenLayersSigma);

    //if(debug && oddLayersMean)
    cout<<"[DTT0CalibrationRMS] Relative T0 sigma for  odd layers "<<oddLayersSigma<<"  even layers"<<evenLayersSigma<<endl;

    //Recompute the mean for odd and even superlayers discarding fluctations
    double oddLayersFinalMean=0;
    double evenLayersFinalMean=0;
    for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
          wiret0 != theRelativeT0PerWire.end();
          ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
          if(((*wiret0).first.layerId().layer()) % 2){
         if(abs((*wiret0).second - oddLayersMean) < (2*oddLayersSigma))
            oddLayersFinalMean = oddLayersFinalMean + (*wiret0).second;
          }
          else{
         if(abs((*wiret0).second - evenLayersMean) < (2*evenLayersSigma))
            evenLayersFinalMean = evenLayersFinalMean + (*wiret0).second;
          }
       }
    }
    oddLayersFinalMean = oddLayersFinalMean/oddLayersDen;
    evenLayersFinalMean = evenLayersFinalMean/evenLayersDen;
    //if(debug && oddLayersMean)
    cout<<"[DTT0CalibrationRMS] Final relative T0 mean for  odd layers "<<oddLayersFinalMean<<"  even layers"<<evenLayersFinalMean<<endl;

    for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
          wiret0 != theRelativeT0PerWire.end();
          ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
          double t0=-999;
          if(((*wiret0).first.layerId().layer()) % 2)
         t0 = (*wiret0).second + (evenLayersFinalMean - oddLayersFinalMean);
          else
         t0 = (*wiret0).second;

          cout << "[DTT0CalibrationRMS] Wire " << (*wiret0).first << " has t0 " << (*wiret0).second
                   << " (relative, after even-odd layer corrections)  "
           << "    sigma " << theSigmaT0PerWire[(*wiret0).first] << endl;

          //Store the results into DB
          t0sRelative->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first],DTTimeUnits::counts); 
       }
    }
     }

     //if(debug)
     cout << "[DTT0CalibrationRMS]Computing relative t0 wrt to chamber average" << endl;
     //Compute the reference for each chamber
     map<DTChamberId,double> sumT0ByChamber;
     map<DTChamberId,int> countT0ByChamber;
     for(DTT0::const_iterator tzero = t0sRelative->begin();
       tzero != t0sRelative->end(); ++tzero) {
        int channelId = tzero->channelId;
        if ( channelId == 0 ) continue;
        DTWireId wireId(channelId);
        DTChamberId chamberId(wireId.chamberId());
        //sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + tzero->t0mean;
        // @@@ better DTT0 usage
        float t0mean_f;
        float t0rms_f;
        t0sRelative->get(wireId,t0mean_f,t0rms_f,DTTimeUnits::counts);
        sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + t0mean_f;
        // @@@ NEW DTT0 END
    countT0ByChamber[chamberId]++;
     }

     //Change reference for each wire and store the new t0s in the new map
     for(DTT0::const_iterator tzero = t0sRelative->begin();
       tzero != t0sRelative->end(); ++tzero) {
    int channelId = tzero->channelId;
    if ( channelId == 0 ) continue;
    DTWireId wireId(channelId);
    DTChamberId chamberId(wireId.chamberId());
    //double t0mean = (tzero->t0mean) - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
    //double t0rms = tzero->t0rms;
    // @@@ better DTT0 usage
    float t0mean_f;
    float t0rms_f;
    t0sRelative->get(wireId,t0mean_f,t0rms_f,DTTimeUnits::counts);
    double t0mean = t0mean_f - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
    double t0rms = t0rms_f;
    // @@@ NEW DTT0 END
    t0sWRTChamber->set(wireId,
          t0mean,
          t0rms,
          DTTimeUnits::counts);
    //if(debug)
    //cout<<"Chamber "<<chamberId<<" has reference "<<(sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
    cout << "Changing t0 of wire " << wireId << " from " << t0mean_f
         << " to " << t0mean << endl;
     }
  }

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

string DTT0CalibrationRMS::getHistoName(const DTWireId& wId) const {
  string histoName = "Ch_" + std::to_string(wId.wheel()) + "_" + std::to_string(wId.station())
                     + "_" + std::to_string(wId.sector()) + "_SL" + std::to_string(wId.superlayer())
                     + "_L" + std::to_string(wId.layer()) + "_W" + std::to_string(wId.wire()) + "_hT0Histo";
  return histoName;
}

string DTT0CalibrationRMS::getHistoName(const DTLayerId& lId) const {
  string histoName = "Ch_" + std::to_string(lId.wheel()) + "_" + std::to_string(lId.station())
                     + "_" + std::to_string(lId.sector()) + "_SL" + std::to_string(lId.superlayer())
                     + "_L" + std::to_string(lId.layer()) + "_hT0Histo";
  return histoName;
}