CTPPS/src/WriteCTPPSPixGainCalibrations.cc

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/ServiceRegistry/interface/Service.h"
0009 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0010 //CTPPS Gain Calibration Conditions Object
0011 #include "CondFormats/PPSObjects/interface/CTPPSPixelGainCalibrations.h"
0012 #include "CondFormats/PPSObjects/interface/CTPPSPixelGainCalibration.h"
0013 #include "CondFormats/PPSObjects/interface/CTPPSPixelIndices.h"
0014 //CTPPS tracker DetId
0015 #include "DataFormats/CTPPSDetId/interface/CTPPSPixelDetId.h"
0016 #include "TFile.h"
0017 #include "TH2F.h"
0018 #include <iostream>
0019 #include <vector>
0020 #include <string>
0021
0022 //
0023 // class declaration
0024 //
0025
0026 class WriteCTPPSPixGainCalibrations : public edm::one::EDAnalyzer<> {
0027 public:
0028   explicit WriteCTPPSPixGainCalibrations(const edm::ParameterSet&);
0029   ~WriteCTPPSPixGainCalibrations() override;
0030
0031   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0032
0033 private:
0034   void beginJob() override;
0035   void analyze(const edm::Event&, const edm::EventSetup&) override;
0036   void endJob() override;
0037   void getHistos();
0038   void fillDB();
0039   void getGainsPedsFromHistos(uint32_t detid,
0040                               int rocId,
0041                               int index,
0042                               std::vector<double>& peds,
0043                               std::vector<double>& gains,
0044                               std::map<int, int>& myindxmap,
0045                               int nrocs);
0046   void setDummyFullPlane(std::vector<float>& peds, std::vector<float>& gains, int npixplane);
0047   // ----------Member data ---------------------------
0048   std::string m_record;
0049   std::string m_inputHistosFileName;
0050   bool m_usedummy;
0051   int npfitmin;
0052   double gainlow, gainhigh;
0053   TFile* m_inputRootFile;
0054   std::map<uint32_t, std::vector<std::string> > detidHistoNameMap;
0055   //  std::map<uint32_t, std::vector<std::string> > detidSlopeNameMap;
0056   //  std::map<uint32_t, std::vector<std::string> > detidInterceptNameMap;
0057   //  std::map<uint32_t, std::vector<std::string> > detidChi2NameMap;
0058   //  std::map<uint32_t, std::vector<std::string> > detidNpfitNameMap;
0059   std::map<uint32_t, std::vector<int> > detidROCsPresent;
0060 };
0061
0062 //
0063 // constructors and destructor
0064 //
0065 WriteCTPPSPixGainCalibrations::WriteCTPPSPixGainCalibrations(const edm::ParameterSet& iConfig)
0066     : m_record(iConfig.getUntrackedParameter<std::string>("record", "CTPPSPixelGainCalibrationsRcd")),
0067       m_inputHistosFileName(iConfig.getUntrackedParameter<std::string>("inputrootfile", "inputfile.root")),
0068       m_usedummy(iConfig.getUntrackedParameter<bool>("useDummyValues", true)),
0069       npfitmin(iConfig.getUntrackedParameter<int>("minimumNpfit", 3)),
0070       gainlow(iConfig.getUntrackedParameter<double>("gainLowLimit", 0.0)),
0071       gainhigh(iConfig.getUntrackedParameter<double>("gainHighLimit", 100.0)) {}
0072
0073 WriteCTPPSPixGainCalibrations::~WriteCTPPSPixGainCalibrations() {
0074   // do anything here that needs to be done at desctruction time
0075   // (e.g. close files, deallocate resources etc.)
0076 }
0077
0078 //
0079 // member functions
0080 //
0081
0082 // ------------ method called for each event  ------------
0083 void WriteCTPPSPixGainCalibrations::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0084   //   using namespace edm;
0085 }
0086
0087 // ------------ method called once each job just before starting event loop  ------------
0088 void WriteCTPPSPixGainCalibrations::beginJob() {}
0089
0090 // ------------ method called once each job just after ending the event loop  ------------
0091 void WriteCTPPSPixGainCalibrations::endJob() {
0092   getHistos();
0093   fillDB();
0094 }
0095
0096 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0097 void WriteCTPPSPixGainCalibrations::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0098   //The following says we do not know what parameters are allowed so do no validation
0099   // Please change this to state exactly what you do use, even if it is no parameters
0100   edm::ParameterSetDescription desc;
0101   desc.setUnknown();
0102   descriptions.addDefault(desc);
0103 }
0104
0105 void WriteCTPPSPixGainCalibrations::getHistos() {
0106   //  std::cout <<"Parsing file " <<m_inputHistosFileName << std::endl;
0107   m_inputRootFile = new TFile(m_inputHistosFileName.c_str());
0108   m_inputRootFile->cd();
0109
0110   int sector[2] = {45, 56};  // arm 0 is sector 45 and arm 1 is sector 56
0111   int nsec = 2;
0112   int station[2] = {0, 2};  // for each arm
0113   int nst = 2;
0114   //int pot[6]={3}; // index of the pot within the 6 pot configuration (vertical top or bottom and horizontal)
0115   int npt = 6;
0116
0117   for (int i = 0; i < nsec; i++)
0118     for (int st = 0; st < nst; st++)
0119       for (int pot = 0; pot < npt; pot++) {
0120         int arm = i;
0121
0122         //Check which pots present
0123         char temppathrp[100];
0124         sprintf(temppathrp, "CTPPS/CTPPS_SEC%d/CTPPS_SEC%d_RP%d%d%d", sector[i], sector[i], arm, station[st], pot);
0125         if (!m_inputRootFile->Get(temppathrp))
0126           continue;
0127
0128         for (int plane = 0; plane < 6; plane++) {
0129           //Check which planes present
0130           char temppathplane[100];
0131           sprintf(temppathplane,
0132                   "CTPPS/CTPPS_SEC%d/CTPPS_SEC%d_RP%d%d%d/CTPPS_SEC%d_RP%d%d%d_PLN%d",
0133                   sector[i],
0134                   sector[i],
0135                   arm,
0136                   station[st],
0137                   pot,
0138                   sector[i],
0139                   arm,
0140                   station[st],
0141                   pot,
0142                   plane);
0143
0144           // do not skip the missing plane, instead put dummy values
0145           //      if(!m_inputRootFile->Get(temppathplane)) continue;
0146
0147           CTPPSPixelDetId mytempid(arm, station[st], pot, plane);
0148           std::vector<std::string> histnamevec;
0149           std::vector<int> listrocs;
0150           for (int roc = 0; roc < 6; roc++) {
0151             char temppathhistos[200];
0152
0153             sprintf(
0154                 temppathhistos,
0155                 "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",
0156                 sector[i],
0157                 sector[i],
0158                 arm,
0159                 station[st],
0160                 pot,
0161                 sector[i],
0162                 arm,
0163                 station[st],
0164                 pot,
0165                 plane,
0166                 sector[i],
0167                 arm,
0168                 station[st],
0169                 pot,
0170                 plane,
0171                 roc);
0172
0173             std::string pathhistos(temppathhistos);
0174             std::string pathslope = pathhistos + "_Slope2D";
0175             std::string pathintercept = pathhistos + "_Intercept2D";
0176             if (m_inputRootFile->Get(pathslope.c_str()) && m_inputRootFile->Get(pathintercept.c_str())) {
0177               histnamevec.push_back(pathhistos);
0178               listrocs.push_back(roc);
0179             }
0180           }
0181           detidHistoNameMap[mytempid.rawId()] = histnamevec;
0182           detidROCsPresent[mytempid.rawId()] = listrocs;
0183           edm::LogInfo("CTPPSPixGainsCalibrationWriter")
0184               << "Raw DetId = " << mytempid.rawId() << " Arm = " << arm << " Sector = " << sector[arm]
0185               << " Station = " << station[st] << " Pot = " << pot << " Plane = " << plane;
0186         }
0187       }
0188 }
0189
0190 void WriteCTPPSPixGainCalibrations::fillDB() {
0191   CTPPSPixelGainCalibrations gainCalibsTest;
0192   CTPPSPixelGainCalibrations gainCalibsTest1;
0193
0194   //  std::cout<<"Here! "<<std::endl;
0195
0196   for (std::map<uint32_t, std::vector<int> >::iterator it = detidROCsPresent.begin(); it != detidROCsPresent.end();
0197        ++it) {
0198     uint32_t tempdetid = it->first;
0199     std::vector<int> rocs = it->second;
0200     unsigned int nrocs = rocs.size();
0201     std::map<int, int> mapIPixIndx;
0202
0203     std::vector<double> gainsFromHistos;
0204     std::vector<double> pedsFromHistos;
0205
0206     CTPPSPixelGainCalibration tempPGCalib;
0207
0208     for (unsigned int i = 0; i < nrocs; i++) {
0209       getGainsPedsFromHistos(tempdetid, i, rocs[i], pedsFromHistos, gainsFromHistos, mapIPixIndx, nrocs);
0210     }
0211
0212     std::vector<float> orderedGains;
0213     std::vector<float> orderedPeds;
0214     for (unsigned int k = 0; k < nrocs * 52 * 80; k++) {
0215       int indx = mapIPixIndx[k];
0216       float tmpped = pedsFromHistos[indx];
0217       float tmpgain = gainsFromHistos[indx];
0218       orderedGains.push_back(tmpgain);
0219       orderedPeds.push_back(tmpped);
0220       tempPGCalib.putData(k, tmpped, tmpgain);
0221     }
0222
0223     if (nrocs == 0) {
0224       edm::LogWarning("CTPPSPixGainsCalibrationWriter") << " plane with detID =" << tempdetid << " is empty";
0225       setDummyFullPlane(orderedPeds, orderedGains, 6 * 52 * 80);
0226     }
0227
0228     gainCalibsTest.setGainCalibration(tempdetid, orderedPeds, orderedGains);
0229     //   std::cout << "Here detid = "<<tempdetid <<std::endl;
0230     gainCalibsTest1.setGainCalibration(tempdetid, tempPGCalib);
0231     //   std::cout << "Here again"<<std::endl;
0232   }
0233   //  std::cout<<" Here 3!"<<std::endl;
0234   edm::Service<cond::service::PoolDBOutputService> mydbservice;
0235   if (!mydbservice.isAvailable()) {
0236     edm::LogError("CTPPSPixGainsCalibrationWriter") << "Db Service Unavailable";
0237     return;
0238   }
0239   mydbservice->writeOneIOV(gainCalibsTest, mydbservice->currentTime(), m_record);
0240 }
0241
0242 void WriteCTPPSPixGainCalibrations::setDummyFullPlane(std::vector<float>& peds,
0243                                                       std::vector<float>& gains,
0244                                                       int npixplane) {
0245   for (int i = 0; i < npixplane; ++i) {
0246     peds.push_back(20.);
0247     gains.push_back(0.5);
0248   }
0249 }
0250
0251 void WriteCTPPSPixGainCalibrations::getGainsPedsFromHistos(uint32_t detid,
0252                                                            int ROCindex,
0253                                                            int rocId,
0254                                                            std::vector<double>& peds,
0255                                                            std::vector<double>& gains,
0256                                                            std::map<int, int>& mymap,
0257                                                            int nrocs) {
0258   CTPPSPixelIndices modulepixels(52 * nrocs / 2, 160);
0259
0260   std::string tmpslopename = detidHistoNameMap[detid][ROCindex] + "_Slope2D";
0261   std::string tmpitcpname = detidHistoNameMap[detid][ROCindex] + "_Intercept2D";
0262   std::string tmpchi2name = detidHistoNameMap[detid][ROCindex] + "_Chisquare2D";
0263   std::string tmpnpfitsname = detidHistoNameMap[detid][ROCindex] + "_Npfits2D";
0264   TH2D* tempslope = (TH2D*)m_inputRootFile->Get(tmpslopename.c_str());
0265   TH2D* tempintrcpt = (TH2D*)m_inputRootFile->Get(tmpitcpname.c_str());
0266   // TH2D * tempchi2    = (TH2D*) m_inputRootFile->Get(tmpchi2name.c_str());
0267   TH2D* tempnpfit = (TH2D*)m_inputRootFile->Get(tmpnpfitsname.c_str());
0268   int ncols = tempslope->GetNbinsX();
0269   int nrows = tempslope->GetNbinsY();
0270   if (nrows != 80 || ncols != 52)
0271     edm::LogWarning("CTPPSPixGainsCalibrationWriter")
0272         << "Something wrong ncols = " << ncols << " and nrows = " << nrows;
0273
0274   for (int jrow = 0; jrow < nrows;
0275        ++jrow)  // when scanning through the 2d histo make sure to avoid underflow bin i or j ==0
0276     for (int icol = 0; icol < ncols; ++icol) {
0277       double tmpslp = tempslope->GetBinContent(icol + 1, jrow + 1);
0278       double tmpgain = (tmpslp == 0.0) ? 0.0 : 1.0 / tmpslp;
0279       double tmpped = tempintrcpt->GetBinContent(icol + 1, jrow + 1);
0280       // check for noisy/badly calibrated pixels
0281       int tmpnpfit = tempnpfit->GetBinContent(icol + 1, jrow + 1);
0282       //double tmpchi2  = tempchi2 -> GetBinContent(icol+1,jrow+1);
0283       if (tmpnpfit < npfitmin || tmpgain < gainlow || tmpgain > gainhigh) {
0284         //  std::cout << " *** Badly calibrated pixel, NPoints = "<<tmpnpfit << " setting dummy values gain = 0.5 and  ped =20.0 ***" <<std::endl;
0285         //  std::cout << " **** bad Pixel column icol = "<<icol <<" and jrow = "<<jrow <<" Name= "<< tmpslopename <<std::endl;
0286         if (m_usedummy) {
0287           tmpgain = 1.0 / 2.0;
0288           tmpped = 20.0;
0289         }
0290         // else  leave as is and set noisy in mask?
0291       }
0292
0293       gains.push_back(tmpgain);
0294       peds.push_back(tmpped);
0295       int modCol = -1;
0296       int modRow = -1;
0297       modulepixels.transformToModule(icol, jrow, rocId, modCol, modRow);
0298       int indx = gains.size() - 1;
0299       int pixIndx = modCol + modRow * (52 * nrocs / 2);
0300       mymap[pixIndx] = indx;
0301     }
0302 }
0303
0304 //define this as a plug-in
0305 DEFINE_FWK_MODULE(WriteCTPPSPixGainCalibrations);