SiStripChannelGainFromDBMiscalibrator.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    CondTools/SiStrip
// Class:      SiStripChannelGainFromDBMiscalibrator
// 
//
// Original Author:  Marco Musich
//         Created:  Tue, 03 Oct 2017 12:57:34 GMT
//
//


// system include files
#include <memory>
#include <iostream>

// user include files
#include "CLHEP/Random/RandGauss.h"
#include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
#include "CalibTracker/Records/interface/SiStripGainRcd.h"
#include "CommonTools/TrackerMap/interface/TrackerMap.h"
#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
#include "CondFormats/DataRecord/interface/SiStripApvGainRcd.h"
#include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
#include "CondFormats/SiStripObjects/interface/SiStripSummary.h"
#include "CondTools/SiStrip/interface/SiStripMiscalibrateHelper.h"
#include "DataFormats/SiStripDetId/interface/SiStripDetId.h" 
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" 
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" 
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h" 
//
// class declaration
//

class SiStripChannelGainFromDBMiscalibrator : public edm::one::EDAnalyzer<>  {
   public:
      explicit SiStripChannelGainFromDBMiscalibrator(const edm::ParameterSet&);
      ~SiStripChannelGainFromDBMiscalibrator() override;

      static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

   private:
      void beginJob() override;
      void analyze(const edm::Event&, const edm::EventSetup&) override;
      std::unique_ptr<SiStripApvGain> getNewObject(const std::map<std::pair<uint32_t,int>,float>& theMap);
      void endJob() override;

      // ----------member data ---------------------------
      const std::string m_Record;  
      const uint32_t m_gainType;
      const bool m_saveMaps;
      const std::vector<edm::ParameterSet> m_parameters;

      std::unique_ptr<TrackerMap> scale_map;
      std::unique_ptr<TrackerMap> smear_map;
      std::unique_ptr<TrackerMap> ratio_map;
      std::unique_ptr<TrackerMap> old_payload_map;
      std::unique_ptr<TrackerMap> new_payload_map;

};

//
// constructors and destructor
//
SiStripChannelGainFromDBMiscalibrator::SiStripChannelGainFromDBMiscalibrator(const edm::ParameterSet& iConfig):
  m_Record{iConfig.getUntrackedParameter<std::string> ("record" , "SiStripApvGainRcd")},
  m_gainType{iConfig.getUntrackedParameter<uint32_t>("gainType",1)},
  m_saveMaps{iConfig.getUntrackedParameter<bool>("saveMaps",true)},
  m_parameters{iConfig.getParameter<std::vector<edm::ParameterSet> >("params")}
{
   //now do what ever initialization is needed
  
  std::string ss_gain = (m_gainType > 0) ? "G2" : "G1";

  scale_map = std::unique_ptr<TrackerMap>(new TrackerMap("scale"));
  scale_map->setTitle("Scale factor averaged by module");
  scale_map->setPalette(1);

  smear_map =std::unique_ptr<TrackerMap>(new TrackerMap("smear"));
  smear_map->setTitle("Smear factor averaged by module");
  smear_map->setPalette(1);

  ratio_map = std::unique_ptr<TrackerMap>(new TrackerMap("ratio"));
  ratio_map->setTitle("Average by module of the "+ss_gain+" Gain payload ratio (new/old)");
  ratio_map->setPalette(1);

  new_payload_map =std::unique_ptr<TrackerMap>(new TrackerMap("new_payload"));
  new_payload_map->setTitle("Tracker Map of Modified "+ss_gain+" Gain payload averaged by module");
  new_payload_map->setPalette(1);
  
  old_payload_map =std::unique_ptr<TrackerMap>(new TrackerMap("old_payload"));
  old_payload_map->setTitle("Tracker Map of Starting "+ss_gain+" Gain Payload averaged by module");
  old_payload_map->setPalette(1);

}


SiStripChannelGainFromDBMiscalibrator::~SiStripChannelGainFromDBMiscalibrator()
{
}


//
// member functions
//

// ------------ method called for each event  ------------
void
SiStripChannelGainFromDBMiscalibrator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;

   edm::ESHandle<TrackerTopology> tTopoHandle;
   iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
   const auto* const tTopo = tTopoHandle.product();

   std::vector<std::string> partitions;

   // fill the list of partitions 
   for(auto& thePSet : m_parameters){
     const std::string partition(thePSet.getParameter<std::string>("partition"));
     // only if it is not yet in the list
     if(std::find(partitions.begin(), partitions.end(), partition) == partitions.end()) {
       partitions.push_back(partition);
     }
   }

   std::map<sistripsummary::TrackerRegion,SiStripMiscalibrate::Smearings> mapOfSmearings;

   for(auto& thePSet : m_parameters){
     
     const std::string partition(thePSet.getParameter<std::string>("partition"));
     sistripsummary::TrackerRegion region = SiStripMiscalibrate::getRegionFromString(partition);
     
     bool    m_doScale(thePSet.getParameter<bool>("doScale"));
     bool    m_doSmear(thePSet.getParameter<bool>("doSmear"));
     double  m_scaleFactor(thePSet.getParameter<double>("scaleFactor"));
     double  m_smearFactor(thePSet.getParameter<double>("smearFactor"));
     
     SiStripMiscalibrate::Smearings params = SiStripMiscalibrate::Smearings();
     params.setSmearing(m_doScale,m_doSmear,m_scaleFactor,m_smearFactor);
     mapOfSmearings[region]=params;
   }
   

   edm::ESHandle<SiStripGain> SiStripApvGain_;
   iSetup.get<SiStripGainRcd>().get(SiStripApvGain_);

   std::map<std::pair<uint32_t,int>,float> theMap,oldPayloadMap;
   
   std::vector<uint32_t> detid;
   SiStripApvGain_->getDetIds(detid);
   for (const auto & d : detid) {
     SiStripApvGain::Range range=SiStripApvGain_->getRange(d,m_gainType);
     float nAPV=0;

     auto regions = SiStripMiscalibrate::getRegionsFromDetId(tTopo,d); 

     // sort by largest to smallest
     std::sort(regions.rbegin(), regions.rend());
     
     SiStripMiscalibrate::Smearings params = SiStripMiscalibrate::Smearings();
     
     for (unsigned int j=0; j<regions.size();j++){
       bool checkRegion = (mapOfSmearings.count(regions[j]) != 0);

       if(!checkRegion) {
     // if the subdetector is not in the list and there's no indication for the whole tracker, just use the default 
     // i.e. no change
     continue;
       } else {
     params = mapOfSmearings[regions[j]];
     break;
       }
     }
     
     scale_map->fill(d,params.m_scaleFactor);
     smear_map->fill(d,params.m_smearFactor);

     for(int it=0;it<range.second-range.first;it++){
       nAPV+=1;
       float Gain=SiStripApvGain_->getApvGain(it,range);
       std::pair<uint32_t,int> index = std::make_pair(d,nAPV);

       oldPayloadMap[index]=Gain;
       
       if(params.m_doScale){
     Gain*=params.m_scaleFactor;
       }
       
       if(params.m_doSmear){
     float smearedGain =  CLHEP::RandGauss::shoot(Gain,params.m_smearFactor);
     Gain=smearedGain;
       }

       theMap[index]=Gain;
       
     } // loop over APVs
   } // loop over DetIds

   std::unique_ptr<SiStripApvGain> theAPVGains = this->getNewObject(theMap);

   // make the payload ratio map
   uint32_t cachedId(0);
   SiStripMiscalibrate::Entry gain_ratio;
   SiStripMiscalibrate::Entry o_gain;
   SiStripMiscalibrate::Entry n_gain;
   for(const auto &element : theMap){

     uint32_t DetId  = element.first.first;
     int      nAPV   = element.first.second;
     float new_gain = element.second;
     float old_gain = oldPayloadMap[std::make_pair(DetId,nAPV)];
     
     // flush the counters
     if(cachedId!=0 && DetId!=cachedId){
       ratio_map->fill(cachedId,gain_ratio.mean());
       old_payload_map->fill(cachedId,o_gain.mean());
       new_payload_map->fill(cachedId,n_gain.mean());

       //auto test = new_payload_map.get()->smoduleMap;

       gain_ratio.reset();
       o_gain.reset();
       n_gain.reset();
     } 
     
     cachedId=DetId;
     gain_ratio.add(new_gain/old_gain);
     o_gain.add(old_gain);
     n_gain.add(new_gain);
   }

   // write out the APVGains record
   edm::Service<cond::service::PoolDBOutputService> poolDbService;
  
   if( poolDbService.isAvailable() )
     poolDbService->writeOne(theAPVGains.get(),poolDbService->currentTime(),m_Record);
   else
     throw std::runtime_error("PoolDBService required.");
 
}

// ------------ method called once each job just before starting event loop  ------------
void 
SiStripChannelGainFromDBMiscalibrator::beginJob()
{
}

// ------------ method called once each job just after ending the event loop  ------------
void 
SiStripChannelGainFromDBMiscalibrator::endJob() 
{
  if(m_saveMaps){

    std::string ss_gain = (m_gainType > 0) ? "G2" : "G1";

    scale_map->save(true,0,0,ss_gain+"_gain_scale_map.pdf");
    scale_map->save(true,0,0,ss_gain+"_gain_scale_map.png"); 
                         
    smear_map->save(true,0,0,ss_gain+"_gain_smear_map.pdf");
    smear_map->save(true,0,0,ss_gain+"_gain_smear_map.png"); 
                         
    ratio_map->save(true,0,0,ss_gain+"_gain_ratio_map.pdf");
    ratio_map->save(true,0,0,ss_gain+"_gain_ratio_map.png"); 

    auto range = SiStripMiscalibrate::getTruncatedRange(old_payload_map.get());

    old_payload_map->save(true,range.first,range.second,"starting_"+ss_gain+"_gain_payload_map.pdf");
    old_payload_map->save(true,range.first,range.second,"starting_"+ss_gain+"_gain_payload_map.png");
                                              
    range = SiStripMiscalibrate::getTruncatedRange(new_payload_map.get());

    new_payload_map->save(true,range.first,range.second,"new_"+ss_gain+"_gain_payload_map.pdf");
    new_payload_map->save(true,range.first,range.second,"new_"+ss_gain+"_gain_payload_map.png");
  }
}

//********************************************************************************//
std::unique_ptr<SiStripApvGain>
SiStripChannelGainFromDBMiscalibrator::getNewObject(const std::map<std::pair<uint32_t,int>,float>& theMap) 
{
  std::unique_ptr<SiStripApvGain> obj = std::unique_ptr<SiStripApvGain>(new SiStripApvGain());
  
  std::vector<float> theSiStripVector;
  uint32_t PreviousDetId = 0; 
  for(const auto &element : theMap){
    uint32_t DetId = element.first.first;
    if(DetId != PreviousDetId){
      if(!theSiStripVector.empty()){
    SiStripApvGain::Range range(theSiStripVector.begin(),theSiStripVector.end());
    if ( !obj->put(PreviousDetId,range) )  printf("Bug to put detId = %i\n",PreviousDetId);
      }
      theSiStripVector.clear();
      PreviousDetId = DetId;
    }
    theSiStripVector.push_back(element.second);
    
    edm::LogInfo("SiStripChannelGainFromDBMiscalibrator")<<" DetId: "<<DetId 
                             <<" APV:   "<<element.first.second
                             <<" Gain:  "<<element.second
                             <<std::endl;
  }
  
  if(!theSiStripVector.empty()){
    SiStripApvGain::Range range(theSiStripVector.begin(),theSiStripVector.end());
    if ( !obj->put(PreviousDetId,range) )  printf("Bug to put detId = %i\n",PreviousDetId);
  }
  
  return obj;
}


// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void
SiStripChannelGainFromDBMiscalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  
  desc.setComment("Creates rescaled / smeared SiStrip Gain payload. Can be used for both G1 and G2."
                  "PoolDBOutputService must be set up for 'SiStripApvGainRcd'.");
  
  edm::ParameterSetDescription descScaler;
  descScaler.setComment("ParameterSet specifying the Strip tracker partition to be scaled / smeared "
                        "by a given factor.");

  descScaler.add<std::string>("partition", "Tracker");
  descScaler.add<bool>("doScale",true);
  descScaler.add<bool>("doSmear",true);
  descScaler.add<double>("scaleFactor", 1.0);
  descScaler.add<double>("smearFactor", 1.0);
  desc.addVPSet("params", descScaler, std::vector<edm::ParameterSet>(1));

  desc.addUntracked<std::string>("record","SiStripApvGainRcd");
  desc.addUntracked<unsigned int>("gainType",1);
  desc.addUntracked<bool>("saveMaps",true);

  descriptions.add("scaleAndSmearSiStripGains", desc);

}

/*--------------------------------------------------------------------*/

//define this as a plug-in
DEFINE_FWK_MODULE(SiStripChannelGainFromDBMiscalibrator);