WriteCTPPSPixGainCalibrations.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
//CTPPS Gain Calibration Conditions Object
#include "CondFormats/PPSObjects/interface/CTPPSPixelGainCalibrations.h"
#include "CondFormats/PPSObjects/interface/CTPPSPixelGainCalibration.h"
#include "CondFormats/PPSObjects/interface/CTPPSPixelIndices.h"
//CTPPS tracker DetId
#include "DataFormats/CTPPSDetId/interface/CTPPSPixelDetId.h"
#include "TFile.h"
#include "TH2F.h"
#include <iostream>
#include <vector>
#include <string>

//
// class declaration
//

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

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

private:
  void beginJob() override;
  void analyze(const edm::Event&, const edm::EventSetup&) override;
  void endJob() override;
  void getHistos();
  void fillDB();
  void getGainsPedsFromHistos(uint32_t detid,
                              int rocId,
                              int index,
                              std::vector<double>& peds,
                              std::vector<double>& gains,
                              std::map<int, int>& myindxmap,
                              int nrocs);
  void setDummyFullPlane(std::vector<float>& peds, std::vector<float>& gains, int npixplane);
  // ----------Member data ---------------------------
  std::string m_record;
  std::string m_inputHistosFileName;
  bool m_usedummy;
  int npfitmin;
  double gainlow, gainhigh;
  TFile* m_inputRootFile;
  std::map<uint32_t, std::vector<std::string> > detidHistoNameMap;
  //  std::map<uint32_t, std::vector<std::string> > detidSlopeNameMap;
  //  std::map<uint32_t, std::vector<std::string> > detidInterceptNameMap;
  //  std::map<uint32_t, std::vector<std::string> > detidChi2NameMap;
  //  std::map<uint32_t, std::vector<std::string> > detidNpfitNameMap;
  std::map<uint32_t, std::vector<int> > detidROCsPresent;
};

//
// constructors and destructor
//
WriteCTPPSPixGainCalibrations::WriteCTPPSPixGainCalibrations(const edm::ParameterSet& iConfig)
    : m_record(iConfig.getUntrackedParameter<std::string>("record", "CTPPSPixelGainCalibrationsRcd")),
      m_inputHistosFileName(iConfig.getUntrackedParameter<std::string>("inputrootfile", "inputfile.root")),
      m_usedummy(iConfig.getUntrackedParameter<bool>("useDummyValues", true)),
      npfitmin(iConfig.getUntrackedParameter<int>("minimumNpfit", 3)),
      gainlow(iConfig.getUntrackedParameter<double>("gainLowLimit", 0.0)),
      gainhigh(iConfig.getUntrackedParameter<double>("gainHighLimit", 100.0)) {}

WriteCTPPSPixGainCalibrations::~WriteCTPPSPixGainCalibrations() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

//
// member functions
//

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

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

// ------------ method called once each job just after ending the event loop  ------------
void WriteCTPPSPixGainCalibrations::endJob() {
  getHistos();
  fillDB();
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void WriteCTPPSPixGainCalibrations::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.setUnknown();
  descriptions.addDefault(desc);
}

void WriteCTPPSPixGainCalibrations::getHistos() {
  //  std::cout <<"Parsing file " <<m_inputHistosFileName << std::endl;
  m_inputRootFile = new TFile(m_inputHistosFileName.c_str());
  m_inputRootFile->cd();

  int sector[2] = {45, 56};  // arm 0 is sector 45 and arm 1 is sector 56
  int nsec = 2;
  int station[2] = {0, 2};  // for each arm
  int nst = 2;
  //int pot[6]={3}; // index of the pot within the 6 pot configuration (vertical top or bottom and horizontal)
  int npt = 6;

  for (int i = 0; i < nsec; i++)
    for (int st = 0; st < nst; st++)
      for (int pot = 0; pot < npt; pot++) {
        int arm = i;

        //Check which pots present
        char temppathrp[100];
        sprintf(temppathrp, "CTPPS/CTPPS_SEC%d/CTPPS_SEC%d_RP%d%d%d", sector[i], sector[i], arm, station[st], pot);
        if (!m_inputRootFile->Get(temppathrp))
          continue;

        for (int plane = 0; plane < 6; plane++) {
          //Check which planes present
          char temppathplane[100];
          sprintf(temppathplane,
                  "CTPPS/CTPPS_SEC%d/CTPPS_SEC%d_RP%d%d%d/CTPPS_SEC%d_RP%d%d%d_PLN%d",
                  sector[i],
                  sector[i],
                  arm,
                  station[st],
                  pot,
                  sector[i],
                  arm,
                  station[st],
                  pot,
                  plane);

          // do not skip the missing plane, instead put dummy values
          //      if(!m_inputRootFile->Get(temppathplane)) continue;

          CTPPSPixelDetId mytempid(arm, station[st], pot, plane);
          std::vector<std::string> histnamevec;
          std::vector<int> listrocs;
          for (int roc = 0; roc < 6; roc++) {
            char temppathhistos[200];

            sprintf(
                temppathhistos,
                "CTPPS/CTPPS_SEC%d/CTPPS_SEC%d_RP%d%d%d/CTPPS_SEC%d_RP%d%d%d_PLN%d/CTPPS_SEC%d_RP%d%d%d_PLN%d_ROC%d",
                sector[i],
                sector[i],
                arm,
                station[st],
                pot,
                sector[i],
                arm,
                station[st],
                pot,
                plane,
                sector[i],
                arm,
                station[st],
                pot,
                plane,
                roc);

            std::string pathhistos(temppathhistos);
            std::string pathslope = pathhistos + "_Slope2D";
            std::string pathintercept = pathhistos + "_Intercept2D";
            if (m_inputRootFile->Get(pathslope.c_str()) && m_inputRootFile->Get(pathintercept.c_str())) {
              histnamevec.push_back(pathhistos);
              listrocs.push_back(roc);
            }
          }
          detidHistoNameMap[mytempid.rawId()] = histnamevec;
          detidROCsPresent[mytempid.rawId()] = listrocs;
          edm::LogInfo("CTPPSPixGainsCalibrationWriter")
              << "Raw DetId = " << mytempid.rawId() << " Arm = " << arm << " Sector = " << sector[arm]
              << " Station = " << station[st] << " Pot = " << pot << " Plane = " << plane;
        }
      }
}

void WriteCTPPSPixGainCalibrations::fillDB() {
  CTPPSPixelGainCalibrations* gainCalibsTest = new CTPPSPixelGainCalibrations();
  CTPPSPixelGainCalibrations* gainCalibsTest1 = new CTPPSPixelGainCalibrations();

  //  std::cout<<"Here! "<<std::endl;

  for (std::map<uint32_t, std::vector<int> >::iterator it = detidROCsPresent.begin(); it != detidROCsPresent.end();
       ++it) {
    uint32_t tempdetid = it->first;
    std::vector<int> rocs = it->second;
    unsigned int nrocs = rocs.size();
    std::map<int, int> mapIPixIndx;

    std::vector<double> gainsFromHistos;
    std::vector<double> pedsFromHistos;

    CTPPSPixelGainCalibration tempPGCalib;

    for (unsigned int i = 0; i < nrocs; i++) {
      getGainsPedsFromHistos(tempdetid, i, rocs[i], pedsFromHistos, gainsFromHistos, mapIPixIndx, nrocs);
    }

    std::vector<float> orderedGains;
    std::vector<float> orderedPeds;
    for (unsigned int k = 0; k < nrocs * 52 * 80; k++) {
      int indx = mapIPixIndx[k];
      float tmpped = pedsFromHistos[indx];
      float tmpgain = gainsFromHistos[indx];
      orderedGains.push_back(tmpgain);
      orderedPeds.push_back(tmpped);
      tempPGCalib.putData(k, tmpped, tmpgain);
    }

    if (nrocs == 0) {
      edm::LogWarning("CTPPSPixGainsCalibrationWriter") << " plane with detID =" << tempdetid << " is empty";
      setDummyFullPlane(orderedPeds, orderedGains, 6 * 52 * 80);
    }

    gainCalibsTest->setGainCalibration(tempdetid, orderedPeds, orderedGains);
    //   std::cout << "Here detid = "<<tempdetid <<std::endl;
    gainCalibsTest1->setGainCalibration(tempdetid, tempPGCalib);
    //   std::cout << "Here again"<<std::endl;
  }
  //  std::cout<<" Here 3!"<<std::endl;
  edm::Service<cond::service::PoolDBOutputService> mydbservice;
  if (!mydbservice.isAvailable()) {
    edm::LogError("CTPPSPixGainsCalibrationWriter") << "Db Service Unavailable";
    return;
  }
  mydbservice->writeOne(gainCalibsTest, mydbservice->currentTime(), m_record);
}

void WriteCTPPSPixGainCalibrations::setDummyFullPlane(std::vector<float>& peds,
                                                      std::vector<float>& gains,
                                                      int npixplane) {
  for (int i = 0; i < npixplane; ++i) {
    peds.push_back(20.);
    gains.push_back(0.5);
  }
}

void WriteCTPPSPixGainCalibrations::getGainsPedsFromHistos(uint32_t detid,
                                                           int ROCindex,
                                                           int rocId,
                                                           std::vector<double>& peds,
                                                           std::vector<double>& gains,
                                                           std::map<int, int>& mymap,
                                                           int nrocs) {
  CTPPSPixelIndices modulepixels(52 * nrocs / 2, 160);

  std::string tmpslopename = detidHistoNameMap[detid][ROCindex] + "_Slope2D";
  std::string tmpitcpname = detidHistoNameMap[detid][ROCindex] + "_Intercept2D";
  std::string tmpchi2name = detidHistoNameMap[detid][ROCindex] + "_Chisquare2D";
  std::string tmpnpfitsname = detidHistoNameMap[detid][ROCindex] + "_Npfits2D";
  TH2D* tempslope = (TH2D*)m_inputRootFile->Get(tmpslopename.c_str());
  TH2D* tempintrcpt = (TH2D*)m_inputRootFile->Get(tmpitcpname.c_str());
  // TH2D * tempchi2    = (TH2D*) m_inputRootFile->Get(tmpchi2name.c_str());
  TH2D* tempnpfit = (TH2D*)m_inputRootFile->Get(tmpnpfitsname.c_str());
  int ncols = tempslope->GetNbinsX();
  int nrows = tempslope->GetNbinsY();
  if (nrows != 80 || ncols != 52)
    edm::LogWarning("CTPPSPixGainsCalibrationWriter")
        << "Something wrong ncols = " << ncols << " and nrows = " << nrows;

  for (int jrow = 0; jrow < nrows;
       ++jrow)  // when scanning through the 2d histo make sure to avoid underflow bin i or j ==0
    for (int icol = 0; icol < ncols; ++icol) {
      double tmpslp = tempslope->GetBinContent(icol + 1, jrow + 1);
      double tmpgain = (tmpslp == 0.0) ? 0.0 : 1.0 / tmpslp;
      double tmpped = tempintrcpt->GetBinContent(icol + 1, jrow + 1);
      // check for noisy/badly calibrated pixels
      int tmpnpfit = tempnpfit->GetBinContent(icol + 1, jrow + 1);
      //double tmpchi2  = tempchi2 -> GetBinContent(icol+1,jrow+1);
      if (tmpnpfit < npfitmin || tmpgain < gainlow || tmpgain > gainhigh) {
        //  std::cout << " *** Badly calibrated pixel, NPoints = "<<tmpnpfit << " setting dummy values gain = 0.5 and  ped =20.0 ***" <<std::endl;
        //  std::cout << " **** bad Pixel column icol = "<<icol <<" and jrow = "<<jrow <<" Name= "<< tmpslopename <<std::endl;
        if (m_usedummy) {
          tmpgain = 1.0 / 2.0;
          tmpped = 20.0;
        }
        // else  leave as is and set noisy in mask?
      }

      gains.push_back(tmpgain);
      peds.push_back(tmpped);
      int modCol = -1;
      int modRow = -1;
      modulepixels.transformToModule(icol, jrow, rocId, modCol, modRow);
      int indx = gains.size() - 1;
      int pixIndx = modCol + modRow * (52 * nrocs / 2);
      mymap[pixIndx] = indx;
    }
}

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