Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:35:42

0001 /**\class EcalPedestalHistory
0002 
0003  Description: <one line class summary>
0004 
0005  Implementation:
0006      <Notes on implementation>
0007 */
0008 //
0009 // $Id: EcalEcalPedestalHistory.cc,v 0.0 2016/05/02 jean fay Exp $
0010 //
0011 //
0012 
0013 #include "EcalPedestalHistory.h"
0014 
0015 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0016 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0017 
0018 //#include<fstream>
0019 
0020 #include "TFile.h"
0021 #include "TTree.h"
0022 #include "TBranch.h"
0023 #include "TH1.h"
0024 #include "TH2.h"
0025 #include "TF1.h"
0026 #include "TProfile.h"
0027 
0028 #include <iostream>
0029 #include <iomanip>
0030 #include <string>
0031 #include <stdexcept>
0032 
0033 using namespace edm;
0034 using namespace cms;
0035 using namespace std;
0036 
0037 //
0038 // constants, enums and typedefs
0039 //
0040 //const Int_t kSample=10;
0041 //
0042 // static data member definitions
0043 //
0044 //int gainValues[3] = {12, 6, 1};
0045 
0046 //
0047 // constructors and destructor
0048 //
0049 
0050 //====================================================================
0051 EcalPedestalHistory::EcalPedestalHistory(const edm::ParameterSet& iConfig) {
0052   //====================================================================
0053   //now do what ever initialization is needed
0054   //  EBDigiCollection_          = iConfig.getParameter<edm::InputTag>("EBDigiCollection");
0055   runnumber_ = iConfig.getUntrackedParameter<int>("runnumber", -1);
0056   ECALType_ = iConfig.getParameter<std::string>("ECALType");
0057   runType_ = iConfig.getParameter<std::string>("runType");
0058   startevent_ = iConfig.getUntrackedParameter<unsigned int>("startevent", 1);
0059 
0060   std::cout << "EcalPedestals Source handler constructor\n" << std::endl;
0061   m_firstRun = static_cast<unsigned int>(atoi(iConfig.getParameter<std::string>("firstRun").c_str()));
0062   m_lastRun = static_cast<unsigned int>(atoi(iConfig.getParameter<std::string>("lastRun").c_str()));
0063   m_sid = iConfig.getParameter<std::string>("OnlineDBSID");
0064   m_user = iConfig.getParameter<std::string>("OnlineDBUser");
0065   m_pass = iConfig.getParameter<std::string>("OnlineDBPassword");
0066   m_locationsource = iConfig.getParameter<std::string>("LocationSource");
0067   m_location = iConfig.getParameter<std::string>("Location");
0068   m_gentag = iConfig.getParameter<std::string>("GenTag");
0069   std::cout << m_sid << "/" << m_user << "/" << m_pass << "/" << m_location << "/" << m_gentag << std::endl;
0070 
0071   vector<int> listDefaults;
0072   listDefaults.push_back(-1);
0073 
0074   cnt_evt_ = 0;
0075   //  cout << "Exiting constructor" << endl;
0076 }  //constructor
0077 
0078 //========================================================================
0079 EcalPedestalHistory::~EcalPedestalHistory() {
0080   //========================================================================
0081   cout << "ANALYSIS FINISHED" << endl;
0082 }  //destructor
0083 
0084 //========================================================================
0085 void EcalPedestalHistory::beginRun(edm::Run const&, edm::EventSetup const& c) {
0086   ///========================================================================
0087 
0088   cout << "Entering beginRun" << endl;
0089   /*     do not use any tag...
0090   edm::ESHandle<EcalChannelStatus> pChannelStatus;
0091   c.get<EcalChannelStatusRcd>().get(pChannelStatus);
0092   const EcalChannelStatus* chStatus = pChannelStatus.product();  
0093   EcalChannelStatusMap::const_iterator chit;
0094   for (int iChannel = 0; iChannel < kEBChannels; iChannel++) {
0095     EBDetId id = EBDetId::unhashIndex(iChannel);
0096     chit = chStatus->getMap().find(id.rawId());
0097     if( chit != chStatus->getMap().end() ) {
0098       EcalChannelStatusCode ch_code = (*chit);
0099       uint16_t statusCode = ch_code.getStatusCode() & 31;
0100       if(statusCode == 1 || (statusCode > 7 && statusCode < 12))
0101     maskedChannels_.push_back(iChannel);
0102     }
0103   }
0104   for (int iChannel = 0; iChannel < kEEChannels; iChannel++) {
0105     EEDetId id = EEDetId::unhashIndex(iChannel);
0106     chit = chStatus->getMap().find(id.rawId());
0107     if( chit != chStatus->getMap().end() ) {
0108       EcalChannelStatusCode ch_code = (*chit);
0109       uint16_t statusCode = ch_code.getStatusCode() & 31;
0110       if(statusCode == 1 || (statusCode > 7 && statusCode < 12))
0111     maskedEEChannels_.push_back(iChannel);
0112     }
0113   }
0114   */
0115   TH1F** hMean = new TH1F*[15];
0116   TH1F** hRMS = new TH1F*[15];
0117   TFile f("PedHist.root", "RECREATE");
0118 
0119   typedef struct {
0120     int iChannel;
0121     int ix;
0122     int iy;
0123     int iz;
0124   } Chan_t;
0125   Chan_t Chan;
0126   Chan.iChannel = -1;
0127   Chan.ix = -1;
0128   Chan.iy = -1;
0129   Chan.iz = -1;
0130 
0131   TTree* tPedChan = new TTree("PedChan", "Channels");  // Output tree for channels
0132   tPedChan->Branch("Channels", &Chan.iChannel, "iChannel/I");
0133   tPedChan->Branch("x", &Chan.ix, "ix/I");
0134   tPedChan->Branch("y", &Chan.iy, "iy/I");
0135   tPedChan->Branch("z", &Chan.iz, "iz/I");
0136   for (int iChannel = 0; iChannel < kEBChannels; iChannel++) {
0137     Chan.iChannel = iChannel;
0138     EBDetId myEBDetId = EBDetId::unhashIndex(iChannel);
0139     Chan.ix = myEBDetId.ieta();  // -85:-1,1:85
0140     Chan.iy = myEBDetId.iphi();  // 1:360
0141     Chan.iz = 0;
0142     if (iChannel % 10000 == 0)
0143       cout << " EB channel " << iChannel << " eta " << Chan.ix << " phi " << Chan.iy << endl;
0144     tPedChan->Fill();
0145   }
0146   for (int iChannel = 0; iChannel < kEEChannels; iChannel++) {
0147     Chan.iChannel = iChannel;
0148     EEDetId myEEDetId = EEDetId::unhashIndex(iChannel);
0149     Chan.ix = myEEDetId.ix();
0150     Chan.iy = myEEDetId.iy();
0151     Chan.iz = myEEDetId.zside();
0152     if (iChannel % 1000 == 0)
0153       cout << " EE channel " << iChannel << " x " << Chan.ix << " y " << Chan.iy << " z " << Chan.iz << endl;
0154     tPedChan->Fill();
0155   }
0156   tPedChan->Write();
0157   tPedChan->Print();
0158 
0159   typedef struct {
0160     int Run;
0161     double Mean[kChannels];
0162     double RMS[kChannels];
0163   } Ped_t;
0164   Ped_t PedVal;
0165   PedVal.Run = -1;  // initialization
0166   for (int iChannel = 0; iChannel < kChannels; iChannel++) {
0167     PedVal.Mean[iChannel] = -1.;
0168     PedVal.RMS[iChannel] = -1.;
0169   }
0170   TTree* tPedHist = new TTree("PedHist", "Pedestal History");  // Output tree for pedestal mean/rms
0171   tPedHist->Branch("Pedestals", &PedVal.Run, "Run/I");
0172   tPedHist->Branch("Mean", PedVal.Mean, "Mean[75848]/D");
0173   tPedHist->Branch("RMS", PedVal.RMS, "RMS[75848]/D");
0174 
0175   // here we retrieve all the runs after the last from online DB
0176   std::cout << "Retrieving run list from ONLINE DB ... " << std::endl;
0177   econn = new EcalCondDBInterface(m_sid, m_user, m_pass);
0178   std::cout << "Connection done" << std::endl;
0179   if (!econn) {
0180     std::cout << " Problem with OMDS: connection parameters " << m_sid << "/" << m_user << "/" << m_pass << std::endl;
0181     throw cms::Exception("OMDS not available");
0182   }
0183 
0184   // these are the online conditions DB classes
0185   RunList my_runlist;
0186   RunTag my_runtag;
0187   LocationDef my_locdef;
0188   RunTypeDef my_rundef;
0189 
0190   my_locdef.setLocation(m_location);
0191   my_rundef.setRunType("PEDESTAL");
0192   my_runtag.setLocationDef(my_locdef);
0193   my_runtag.setRunTypeDef(my_rundef);
0194   my_runtag.setGeneralTag(m_gentag);
0195 
0196   // here we retrieve the Monitoring run records
0197   MonVersionDef monverdef;
0198   monverdef.setMonitoringVersion("test01");
0199   MonRunTag mon_tag;
0200   //    mon_tag.setGeneralTag("CMSSW");
0201   mon_tag.setGeneralTag("CMSSW-offline-private");
0202   mon_tag.setMonVersionDef(monverdef);
0203   MonRunList mon_list;
0204   mon_list.setMonRunTag(mon_tag);
0205   mon_list.setRunTag(my_runtag);
0206   //    mon_list=econn->fetchMonRunList(my_runtag, mon_tag);
0207   unsigned int min_run = 0, max_since = 0;
0208   if (m_firstRun < max_since) {
0209     min_run = max_since + 1;  // we have to add 1 to the last transferred one
0210   } else {
0211     min_run = m_firstRun;
0212   }
0213 
0214   unsigned int max_run = m_lastRun;
0215   mon_list = econn->fetchMonRunList(my_runtag, mon_tag, min_run, max_run);
0216 
0217   std::vector<MonRunIOV> mon_run_vec = mon_list.getRuns();
0218   int mon_runs = mon_run_vec.size();
0219   std::cout << "number of Mon runs is : " << mon_runs << std::endl;
0220 
0221   if (mon_runs > 0) {
0222     int NbChan = 0;
0223     for (int iChannel = 0; iChannel < kEBChannels; iChannel++) {
0224       if (iChannel % 10000 == 1) {
0225         hMean[NbChan] = new TH1F(Form("Mean_%i", NbChan), Form("Mean EB %i", iChannel), mon_runs, 0., mon_runs);
0226         hRMS[NbChan] = new TH1F(Form("RMS_%i", NbChan), Form("RMS EB %i", iChannel), mon_runs, 0., mon_runs);
0227         NbChan++;
0228       }
0229     }
0230     for (int iChannel = 0; iChannel < kEEChannels; iChannel++) {
0231       if (iChannel % 2000 == 1) {
0232         hMean[NbChan] = new TH1F(Form("Mean_%i", NbChan), Form("Mean EE %i", iChannel), mon_runs, 0., mon_runs);
0233         hRMS[NbChan] = new TH1F(Form("RMS_%i", NbChan), Form("RMS EE %i", iChannel), mon_runs, 0., mon_runs);
0234         NbChan++;
0235       }
0236     }
0237 
0238     //    int krmax = std::min(mon_runs, 30);
0239     int krmax = mon_runs;
0240     for (int kr = 0; kr < krmax; kr++) {
0241       std::cout << "-kr------:  " << kr << std::endl;
0242 
0243       unsigned int irun = static_cast<unsigned int>(mon_run_vec[kr].getRunIOV().getRunNumber());
0244       std::cout << "retrieve the data for run number: " << irun << std::endl;
0245       if (mon_run_vec[kr].getSubRunNumber() <= 1) {
0246         // retrieve the data for a given run
0247         RunIOV runiov_prime = mon_run_vec[kr].getRunIOV();
0248         // retrieve the pedestals from OMDS for this run
0249         std::map<EcalLogicID, MonPedestalsDat> dataset_mon;
0250         econn->fetchDataSet(&dataset_mon, &mon_run_vec[kr]);
0251         std::cout << "OMDS record for run " << irun << " is made of " << dataset_mon.size() << std::endl;
0252         int nEB = 0, nEE = 0;
0253         typedef std::map<EcalLogicID, MonPedestalsDat>::const_iterator CImon;
0254         EcalLogicID ecid_xt;
0255         MonPedestalsDat rd_ped;
0256 
0257         // this to validate ...
0258         int nbad = 0;
0259         for (CImon p = dataset_mon.begin(); p != dataset_mon.end(); p++) {
0260           ecid_xt = p->first;
0261           rd_ped = p->second;
0262           int sm_num = ecid_xt.getID1();
0263           int xt_num = ecid_xt.getID2();
0264           int yt_num = ecid_xt.getID3();
0265 
0266           //checkPedestal
0267           bool result = true;
0268           if (rd_ped.getPedRMSG12() > 3 || rd_ped.getPedRMSG12() <= 0 || rd_ped.getPedRMSG6() > 2 ||
0269               rd_ped.getPedRMSG12() <= 0 || rd_ped.getPedRMSG1() > 1 || rd_ped.getPedRMSG1() <= 0 ||
0270               rd_ped.getPedMeanG12() > 300 || rd_ped.getPedMeanG12() <= 100 || rd_ped.getPedMeanG6() > 300 ||
0271               rd_ped.getPedMeanG6() <= 100 || rd_ped.getPedMeanG1() > 300 || rd_ped.getPedMeanG6() <= 100)
0272             result = false;
0273 
0274           // here we check and count how many bad channels we have
0275           if (!result) {
0276             nbad++;
0277             if (nbad < 10)
0278               std::cout << "BAD LIST: channel " << sm_num << "/" << xt_num << "/" << yt_num << "ped/rms "
0279                         << rd_ped.getPedMeanG12() << "/" << rd_ped.getPedRMSG12() << std::endl;
0280           }
0281           if (ecid_xt.getName() == "EB_crystal_number") {
0282             nEB++;
0283           } else {
0284             nEE++;
0285           }
0286         }  // end loop over pedestal data
0287         // ok or bad? A bad run is for more than 5% bad channels
0288 
0289         //        if(nbad<(dataset_mon.size()*0.1)){
0290         if (nbad < (dataset_mon.size() * 0.05) &&
0291             (nEB > 10200 || nEE > 2460)) {  // this is good run, fill histo and tree
0292           PedVal.Run = irun;
0293           int NbChan = 0;
0294           for (CImon p = dataset_mon.begin(); p != dataset_mon.end(); p++) {
0295             ecid_xt = p->first;
0296             rd_ped = p->second;
0297             int sm_num = ecid_xt.getID1();
0298             int xt_num = ecid_xt.getID2();
0299             int yt_num = ecid_xt.getID3();
0300 
0301             if (ecid_xt.getName() == "EB_crystal_number") {  // Barrel
0302               EBDetId ebdetid(sm_num, xt_num, EBDetId::SMCRYSTALMODE);
0303               int iChannel = ebdetid.hashedIndex();
0304               if (iChannel < 0 || iChannel > 61200)
0305                 cout << " SM " << sm_num << " Chan in SM " << xt_num << " IChannel " << iChannel << endl;
0306               if (iChannel % 10000 == 1) {
0307                 hMean[NbChan]->Fill(kr, rd_ped.getPedMeanG12());
0308                 hRMS[NbChan]->Fill(kr, rd_ped.getPedRMSG12());
0309                 NbChan++;
0310               }
0311               PedVal.Mean[iChannel] = rd_ped.getPedMeanG12();
0312               PedVal.RMS[iChannel] = rd_ped.getPedRMSG12();
0313               if (iChannel % 10000 == 0)
0314                 cout << " channel " << iChannel << " mean " << PedVal.Mean[iChannel] << " RMS " << PedVal.RMS[iChannel]
0315                      << endl;
0316             } else {  // Endcaps
0317               if (EEDetId::validDetId(xt_num, yt_num, sm_num)) {
0318                 EEDetId eedetid(xt_num, yt_num, sm_num);
0319                 int iChannel = eedetid.hashedIndex();
0320                 if (iChannel < 0 || iChannel > 14648)
0321                   cout << " x " << sm_num << " y " << xt_num << " z " << yt_num << " IChannel " << iChannel << endl;
0322                 if (iChannel % 2000 == 1) {
0323                   hMean[NbChan]->Fill(kr, rd_ped.getPedMeanG12());
0324                   hRMS[NbChan]->Fill(kr, rd_ped.getPedRMSG12());
0325                   NbChan++;
0326                 }
0327                 int iChanEE = kEBChannels + iChannel;
0328                 //      cout << " channel EE " << iChanEE << endl;
0329                 PedVal.Mean[iChanEE] = rd_ped.getPedMeanG12();
0330                 PedVal.RMS[iChanEE] = rd_ped.getPedRMSG12();
0331               }  // valid ee Id
0332             }  // Endcaps
0333           }  // loop over channels
0334           tPedHist->Fill();
0335           cout << " We got a good run " << irun << endl;
0336         }  // good run
0337       }  // mon_run_vec
0338     }  // loop over all runs
0339   }  // number of runs > 0
0340   cout << "Exiting beginRun" << endl;
0341   for (int NbChan = 0; NbChan < 15; NbChan++) {
0342     if (hMean[NbChan]->GetEntries() > 0.) {  // save only when filled!
0343       hMean[NbChan]->Write();
0344       hRMS[NbChan]->Write();
0345     }
0346   }
0347   tPedHist->Write();
0348   tPedHist->Print();
0349   f.Close();
0350 
0351 }  //beginRun
0352 
0353 //========================================================================
0354 void EcalPedestalHistory::endRun(edm::Run const&, edm::EventSetup const& c) {
0355   //========================================================================
0356 
0357 }  //endRun
0358 
0359 //========================================================================
0360 void EcalPedestalHistory::beginJob() {
0361   ///========================================================================
0362 
0363 }  //beginJob
0364 
0365 //========================================================================
0366 void EcalPedestalHistory::endJob() {
0367   //========================================================================
0368 
0369 }  //endJob
0370 
0371 //
0372 // member functions
0373 //
0374 
0375 //========================================================================
0376 void EcalPedestalHistory::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0377   //========================================================================
0378 
0379   if (cnt_evt_ == 0) {
0380     if (ECALType_ == "EB" || ECALType_ == "EA") {
0381       cout << " Barrel data : nb channels " << kEBChannels << endl;
0382     } else if (ECALType_ == "EE" || ECALType_ == "EA") {
0383       cout << " End cap data : nb channels " << kEEChannels << endl;
0384     } else {
0385       cout << " strange ECALtype : " << ECALType_ << " abort " << endl;
0386       return;
0387     }
0388     /*
0389     int NbOfmaskedChannels =  maskedChannels_.size();
0390     cout << " Nb masked EB channels " << NbOfmaskedChannels << endl;
0391     for (vector<int>::iterator iter = maskedChannels_.begin(); iter != maskedChannels_.end(); ++iter)
0392       cout<< " : masked channel " << *(iter) << endl;
0393     NbOfmaskedChannels =  maskedEEChannels_.size();
0394     cout << " Nb masked EE channels " << NbOfmaskedChannels << endl;
0395     for (vector<int>::iterator iter = maskedEEChannels_.begin(); iter != maskedEEChannels_.end(); ++iter)
0396       cout<< " : masked channel " << *(iter) << endl;
0397     */
0398   }
0399   cnt_evt_++;
0400 
0401 }  //analyze
0402 
0403 //define this as a plug-in
0404 DEFINE_FWK_MODULE(EcalPedestalHistory);