Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:25 

0001 /**_________________________________________________________________
0002 class:   CorrPCCProducer.cc
0003 
0004 description: Computes the type 1 and type 2 corrections to the luminosity
0005                 type 1 - first (spillover from previous BXs real clusters)
0006                 type 2 - after (comes from real activation)
0007                 pedestal - is a constant noise term for low lumi period
0008 
0009 authors:Sam Higginbotham (shigginb@cern.ch) and Chris Palmer (capalmer@cern.ch) , Jose Benitez (jose.benitez@cern.ch)
0010 
0011 ________________________________________________________________**/
0012 #include <memory>
0013 #include <string>
0014 #include <vector>
0015 #include <boost/serialization/vector.hpp>
0016 #include <iostream>
0017 #include <map>
0018 #include <utility>
0019 #include "DQMServices/Core/interface/DQMOneEDAnalyzer.h"
0020 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0021 #include "CondFormats/Luminosity/interface/LumiCorrections.h"
0022 #include "CondFormats/DataRecord/interface/LumiCorrectionsRcd.h"
0023 #include "CondFormats/Serialization/interface/Serializable.h"
0024 #include "DataFormats/Luminosity/interface/PixelClusterCounts.h"
0025 #include "DataFormats/Luminosity/interface/LumiInfo.h"
0026 #include "DataFormats/Luminosity/interface/LumiConstants.h"
0027 #include "DataFormats/Provenance/interface/LuminosityBlockRange.h"
0028 #include "DataFormats/Common/interface/Handle.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/Framework/interface/ConsumesCollector.h"
0032 #include "FWCore/Framework/interface/FileBlock.h"
0033 #include "FWCore/Framework/interface/ESHandle.h"
0034 #include "FWCore/Framework/interface/EventSetup.h"
0035 #include "FWCore/Framework/interface/IOVSyncValue.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/one/EDProducer.h"
0038 #include "FWCore/Framework/interface/Event.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0040 #include "FWCore/Utilities/interface/EDGetToken.h"
0041 #include "FWCore/ServiceRegistry/interface/Service.h"
0042 #include "FWCore/Framework/interface/LuminosityBlock.h"
0043 #include "FWCore/Framework/interface/Run.h"
0044 #include "TMath.h"
0045 #include "TH1.h"
0046 #include "TGraph.h"
0047 #include "TGraphErrors.h"
0048 #include "TFile.h"
0049 
0050 class CorrPCCProducer : public DQMOneEDAnalyzer<edm::one::WatchLuminosityBlocks> {
0051 public:
0052   explicit CorrPCCProducer(const edm::ParameterSet&);
0053   ~CorrPCCProducer() override;
0054 
0055 private:
0056   void beginLuminosityBlock(edm::LuminosityBlock const& lumiSeg, const edm::EventSetup& iSetup) final;
0057   void endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, const edm::EventSetup& iSetup) final;
0058   void dqmEndRun(edm::Run const& runSeg, const edm::EventSetup& iSetup) final;
0059   void dqmEndRunProduce(const edm::Run& runSeg, const edm::EventSetup& iSetup);
0060   void endJob() final;
0061 
0062   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0063 
0064   void makeCorrectionTemplate();
0065   float getMaximum(std::vector<float>);
0066   void estimateType1Frac(std::vector<float>, float&);
0067   void evaluateCorrectionResiduals(std::vector<float>);
0068   void calculateCorrections(std::vector<float>, std::vector<float>&, float&);
0069   void resetBlock();
0070 
0071   edm::EDGetTokenT<LumiInfo> lumiInfoToken;
0072   std::string pccSrc_;    //input file EDproducer module label
0073   std::string prodInst_;  //input file product instance
0074 
0075   std::vector<float> rawlumiBX_;         //new vector containing clusters per bxid
0076   std::vector<float> errOnLumiByBX_;     //standard error per bx
0077   std::vector<float> totalLumiByBX_;     //summed lumi
0078   std::vector<float> totalLumiByBXAvg_;  //summed lumi
0079   std::vector<float> events_;            //Number of events in each BX
0080   std::vector<float> correctionTemplate_;
0081   std::vector<float> correctionScaleFactors_;  //list of scale factors to apply.
0082   float overallCorrection_;                    //The Overall correction to the integrated luminosity
0083 
0084   unsigned int iBlock = 0;
0085   unsigned int minimumNumberOfEvents;
0086 
0087   std::map<unsigned int, LumiInfo*> lumiInfoMapPerLS;
0088   std::vector<unsigned int> lumiSections;
0089   std::map<std::pair<unsigned int, unsigned int>, LumiInfo*>::iterator lumiInfoMapIterator;
0090   std::map<std::pair<unsigned int, unsigned int>, LumiInfo*>
0091       lumiInfoMap;  //map to obtain iov for lumiOb corrections to the luminosity.
0092   std::map<std::pair<unsigned int, unsigned int>, unsigned int> lumiInfoCounter;  // number of lumiSections in this block
0093 
0094   TH1F* corrlumiAvg_h;
0095   TH1F* scaleFactorAvg_h;
0096   TH1F* lumiAvg_h;
0097   TH1F* type1FracHist;
0098   TH1F* type1resHist;
0099   TH1F* type2resHist;
0100 
0101   unsigned int maxLS = 3500;
0102   MonitorElement* Type1FracMon;
0103   MonitorElement* Type1ResMon;
0104   MonitorElement* Type2ResMon;
0105 
0106   TGraphErrors* type1FracGraph;
0107   TGraphErrors* type1resGraph;
0108   TGraphErrors* type2resGraph;
0109   TList* hlist;  //list for the clusters and corrections
0110   TFile* histoFile;
0111 
0112   float type1Frac;
0113   float mean_type1_residual;          //Type 1 residual
0114   float mean_type2_residual;          //Type 2 residual
0115   float mean_type1_residual_unc;      //Type 1 residual uncertainty rms
0116   float mean_type2_residual_unc;      //Type 2 residual uncertainty rms
0117   unsigned int nTrain;                //Number of bunch trains used in calc type 1 and 2 res, frac.
0118   unsigned int countLumi_;            //The lumisection count... the size of the lumiblock
0119   unsigned int approxLumiBlockSize_;  //The number of lumisections per block.
0120   unsigned int thisLS;                //Ending lumisection for the iov that we save with the lumiInfo object.
0121 
0122   double type2_a_;  //amplitude for the type 2 correction
0123   double type2_b_;  //decay width for the type 2 correction
0124 
0125   float pedestal;
0126   float pedestal_unc;
0127   TGraphErrors* pedestalGraph;
0128 
0129   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0130 };
0131 
0132 //--------------------------------------------------------------------------------------------------
0133 CorrPCCProducer::CorrPCCProducer(const edm::ParameterSet& iConfig) {
0134   pccSrc_ =
0135       iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<std::string>("inLumiObLabel");
0136   prodInst_ =
0137       iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<std::string>("ProdInst");
0138   approxLumiBlockSize_ =
0139       iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<int>("approxLumiBlockSize");
0140   type2_a_ = iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<double>("type2_a");
0141   type2_b_ = iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<double>("type2_b");
0142   countLumi_ = 0;
0143   minimumNumberOfEvents = 1000;
0144 
0145   totalLumiByBX_.resize(LumiConstants::numBX);
0146   totalLumiByBXAvg_.resize(LumiConstants::numBX);
0147   events_.resize(LumiConstants::numBX);
0148   correctionScaleFactors_.resize(LumiConstants::numBX);
0149   correctionTemplate_.resize(LumiConstants::numBX);
0150 
0151   resetBlock();
0152 
0153   makeCorrectionTemplate();
0154 
0155   edm::InputTag inputPCCTag_(pccSrc_, prodInst_);
0156 
0157   lumiInfoToken = consumes<LumiInfo, edm::InLumi>(inputPCCTag_);
0158 
0159   histoFile = new TFile("CorrectionHisto.root", "RECREATE");
0160 
0161   type1FracGraph = new TGraphErrors();
0162   type1resGraph = new TGraphErrors();
0163   type2resGraph = new TGraphErrors();
0164   type1FracGraph->SetName("Type1Fraction");
0165   type1resGraph->SetName("Type1Res");
0166   type2resGraph->SetName("Type2Res");
0167   type1FracGraph->GetYaxis()->SetTitle("Type 1 Fraction");
0168   type1resGraph->GetYaxis()->SetTitle("Type 1 Residual");
0169   type2resGraph->GetYaxis()->SetTitle("Type 2 Residual");
0170   type1FracGraph->GetXaxis()->SetTitle("Unique LS ID");
0171   type1resGraph->GetXaxis()->SetTitle("Unique LS ID");
0172   type2resGraph->GetXaxis()->SetTitle("Unique LS ID");
0173   type1FracGraph->SetMarkerStyle(8);
0174   type1resGraph->SetMarkerStyle(8);
0175   type2resGraph->SetMarkerStyle(8);
0176 
0177   pedestalGraph = new TGraphErrors();
0178   pedestalGraph->SetName("Pedestal");
0179   pedestalGraph->GetYaxis()->SetTitle("pedestal value (counts) per lumi section");
0180   pedestalGraph->GetXaxis()->SetTitle("Unique LS ID");
0181   pedestalGraph->SetMarkerStyle(8);
0182 }
0183 
0184 //--------------------------------------------------------------------------------------------------
0185 CorrPCCProducer::~CorrPCCProducer() {}
0186 
0187 //--------------------------------------------------------------------------------------------------
0188 // This method builds the single bunch response given an exponential function (type 2 only)
0189 void CorrPCCProducer::makeCorrectionTemplate() {
0190   for (unsigned int bx = 1; bx < LumiConstants::numBX; bx++) {
0191     correctionTemplate_.at(bx) = type2_a_ * exp(-(float(bx) - 1) * type2_b_);
0192   }
0193 }
0194 
0195 //--------------------------------------------------------------------------------------------------
0196 // Finds max lumi value
0197 float CorrPCCProducer::getMaximum(std::vector<float> lumi_vector) {
0198   float max_lumi = 0;
0199   for (size_t i = 0; i < lumi_vector.size(); i++) {
0200     if (lumi_vector.at(i) > max_lumi)
0201       max_lumi = lumi_vector.at(i);
0202   }
0203   return max_lumi;
0204 }
0205 
0206 //--------------------------------------------------------------------------------------------------
0207 // This method takes luminosity from the last bunch in a train and makes a comparison with
0208 // the follow non-active bunch crossing to estimate the spill over fraction (type 1 afterglow).
0209 void CorrPCCProducer::estimateType1Frac(std::vector<float> uncorrPCCPerBX, float& type1Frac) {
0210   std::vector<float> corrected_tmp_;
0211   for (size_t i = 0; i < uncorrPCCPerBX.size(); i++) {
0212     corrected_tmp_.push_back(uncorrPCCPerBX.at(i));
0213   }
0214 
0215   //Apply initial type 1 correction
0216   for (size_t k = 0; k < LumiConstants::numBX - 1; k++) {
0217     float bin_k = corrected_tmp_.at(k);
0218     corrected_tmp_.at(k + 1) = corrected_tmp_.at(k + 1) - type1Frac * bin_k;
0219   }
0220 
0221   //Apply type 2 correction
0222   for (size_t i = 0; i < LumiConstants::numBX - 1; i++) {
0223     for (size_t j = i + 1; j < i + LumiConstants::numBX - 1; j++) {
0224       float bin_i = corrected_tmp_.at(i);
0225       if (j < LumiConstants::numBX) {
0226         corrected_tmp_.at(j) = corrected_tmp_.at(j) - bin_i * correctionTemplate_.at(j - i);
0227       } else {
0228         corrected_tmp_.at(j - LumiConstants::numBX) =
0229             corrected_tmp_.at(j - LumiConstants::numBX) - bin_i * correctionTemplate_.at(j - i);
0230       }
0231     }
0232   }
0233 
0234   //Apply additional iteration for type 1 correction
0235   evaluateCorrectionResiduals(corrected_tmp_);
0236   type1Frac += mean_type1_residual;
0237 }
0238 
0239 //--------------------------------------------------------------------------------------------------
0240 void CorrPCCProducer::evaluateCorrectionResiduals(std::vector<float> corrected_tmp_) {
0241   float lumiMax = getMaximum(corrected_tmp_);
0242   float threshold = lumiMax * 0.2;
0243 
0244   mean_type1_residual = 0;
0245   mean_type2_residual = 0;
0246   nTrain = 0;
0247   float lumi = 0;
0248   std::vector<float> afterGlow;
0249   TH1F type1("type1", "", 1000, -0.5, 0.5);
0250   TH1F type2("type2", "", 1000, -0.5, 0.5);
0251   for (size_t ibx = 2; ibx < LumiConstants::numBX - 5; ibx++) {
0252     lumi = corrected_tmp_.at(ibx);
0253     afterGlow.clear();
0254     afterGlow.push_back(corrected_tmp_.at(ibx + 1));
0255     afterGlow.push_back(corrected_tmp_.at(ibx + 2));
0256 
0257     //Where type 1 and type 2 residuals are computed
0258     if (lumi > threshold && afterGlow[0] < threshold && afterGlow[1] < threshold) {
0259       for (int index = 3; index < 6; index++) {
0260         float thisAfterGlow = corrected_tmp_.at(ibx + index);
0261         if (thisAfterGlow < threshold) {
0262           afterGlow.push_back(thisAfterGlow);
0263         } else {
0264           break;
0265         }
0266       }
0267       float thisType1 = 0;
0268       float thisType2 = 0;
0269       if (afterGlow.size() > 1) {
0270         int nAfter = 0;
0271         for (unsigned int index = 1; index < afterGlow.size(); index++) {
0272           thisType2 += afterGlow[index];
0273           type2.Fill(afterGlow[index] / lumi);
0274           nAfter++;
0275         }
0276         thisType2 /= nAfter;
0277       }
0278       thisType1 = (afterGlow[0] - thisType2) / lumi;
0279 
0280       type1.Fill(thisType1);
0281       nTrain += 1;
0282     }
0283   }
0284 
0285   mean_type1_residual = type1.GetMean();  //Calculate the mean value of the type 1 residual
0286   mean_type2_residual = type2.GetMean();  //Calculate the mean value of the type 2 residual
0287   mean_type1_residual_unc = type1.GetMeanError();
0288   mean_type2_residual_unc = type2.GetMeanError();
0289 
0290   histoFile->cd();
0291   type1.Write();
0292   type2.Write();
0293 }
0294 
0295 //--------------------------------------------------------------------------------------------------
0296 void CorrPCCProducer::calculateCorrections(std::vector<float> uncorrected,
0297                                            std::vector<float>& correctionScaleFactors_,
0298                                            float& overallCorrection_) {
0299   type1Frac = 0;
0300 
0301   int nTrials = 4;
0302 
0303   for (int trial = 0; trial < nTrials; trial++) {
0304     estimateType1Frac(uncorrected, type1Frac);
0305     edm::LogInfo("INFO") << "type 1 fraction after iteration " << trial << " is  " << type1Frac;
0306   }
0307 
0308   //correction should never be negative
0309   type1Frac = std::max(0.0, (double)type1Frac);
0310 
0311   std::vector<float> corrected_tmp_;
0312   for (size_t i = 0; i < uncorrected.size(); i++) {
0313     corrected_tmp_.push_back(uncorrected.at(i));
0314   }
0315 
0316   //Apply all corrections
0317   for (size_t i = 0; i < LumiConstants::numBX - 1; i++) {
0318     // type 1 - first (spillover from previous BXs real clusters)
0319     float bin_i = corrected_tmp_.at(i);
0320     corrected_tmp_.at(i + 1) = corrected_tmp_.at(i + 1) - type1Frac * bin_i;
0321 
0322     // type 2 - after (comes from real activation)
0323     bin_i = corrected_tmp_.at(i);
0324     for (size_t j = i + 1; j < i + LumiConstants::numBX - 1; j++) {
0325       if (j < LumiConstants::numBX) {
0326         corrected_tmp_.at(j) = corrected_tmp_.at(j) - bin_i * correctionTemplate_.at(j - i);
0327       } else {
0328         corrected_tmp_.at(j - LumiConstants::numBX) =
0329             corrected_tmp_.at(j - LumiConstants::numBX) - bin_i * correctionTemplate_.at(j - i);
0330       }
0331     }
0332   }
0333 
0334   float lumiMax = getMaximum(corrected_tmp_);
0335   float threshold = lumiMax * 0.2;
0336 
0337   //here subtract the pedestal
0338   pedestal = 0.;
0339   pedestal_unc = 0.;
0340   int nped = 0;
0341   for (size_t i = 0; i < LumiConstants::numBX; i++) {
0342     if (corrected_tmp_.at(i) < threshold) {
0343       pedestal += corrected_tmp_.at(i);
0344       nped++;
0345     }
0346   }
0347   if (nped > 0) {
0348     pedestal_unc = sqrt(pedestal) / nped;
0349     pedestal = pedestal / nped;
0350   }
0351   for (size_t i = 0; i < LumiConstants::numBX; i++) {
0352     corrected_tmp_.at(i) = corrected_tmp_.at(i) - pedestal;
0353   }
0354 
0355   evaluateCorrectionResiduals(corrected_tmp_);
0356 
0357   float integral_uncorr_clusters = 0;
0358   float integral_corr_clusters = 0;
0359 
0360   //Calculate Per-BX correction factor and overall correction factor
0361   for (size_t ibx = 0; ibx < corrected_tmp_.size(); ibx++) {
0362     if (corrected_tmp_.at(ibx) > threshold) {
0363       integral_uncorr_clusters += uncorrected.at(ibx);
0364       integral_corr_clusters += corrected_tmp_.at(ibx);
0365     }
0366     if (corrected_tmp_.at(ibx) != 0.0 && uncorrected.at(ibx) != 0.0) {
0367       correctionScaleFactors_.at(ibx) = corrected_tmp_.at(ibx) / uncorrected.at(ibx);
0368     } else {
0369       correctionScaleFactors_.at(ibx) = 0.0;
0370     }
0371   }
0372 
0373   overallCorrection_ = integral_corr_clusters / integral_uncorr_clusters;
0374 }
0375 
0376 //--------------------------------------------------------------------------------------------------
0377 void CorrPCCProducer::beginLuminosityBlock(edm::LuminosityBlock const& lumiSeg, const edm::EventSetup& iSetup) {
0378   countLumi_++;
0379 }
0380 
0381 //--------------------------------------------------------------------------------------------------
0382 void CorrPCCProducer::endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, const edm::EventSetup& iSetup) {
0383   thisLS = lumiSeg.luminosityBlock();
0384 
0385   edm::Handle<LumiInfo> PCCHandle;
0386   lumiSeg.getByToken(lumiInfoToken, PCCHandle);
0387 
0388   const LumiInfo& inLumiOb = *(PCCHandle.product());
0389 
0390   errOnLumiByBX_ = inLumiOb.getErrorLumiAllBX();
0391 
0392   unsigned int totalEvents = 0;
0393   for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
0394     totalEvents += errOnLumiByBX_[bx];
0395   }
0396 
0397   if (totalEvents < minimumNumberOfEvents) {
0398     edm::LogInfo("INFO") << "number of events in this LS is too few " << totalEvents;
0399     //return;
0400   } else {
0401     edm::LogInfo("INFO") << "Skipping Lumisection " << thisLS;
0402   }
0403 
0404   lumiInfoMapPerLS[thisLS] = new LumiInfo();
0405   totalLumiByBX_ = inLumiOb.getInstLumiAllBX();
0406   events_ = inLumiOb.getErrorLumiAllBX();
0407   lumiInfoMapPerLS[thisLS]->setInstLumiAllBX(totalLumiByBX_);
0408   lumiInfoMapPerLS[thisLS]->setErrorLumiAllBX(events_);
0409   lumiSections.push_back(thisLS);
0410 }
0411 
0412 //--------------------------------------------------------------------------------------------------
0413 void CorrPCCProducer::dqmEndRun(edm::Run const& runSeg, const edm::EventSetup& iSetup) {
0414   // TODO: why was this code not put here in the first place?
0415   dqmEndRunProduce(runSeg, iSetup);
0416 }
0417 
0418 //--------------------------------------------------------------------------------------------------
0419 void CorrPCCProducer::dqmEndRunProduce(edm::Run const& runSeg, const edm::EventSetup& iSetup) {
0420   if (lumiSections.empty()) {
0421     return;
0422   }
0423 
0424   std::sort(lumiSections.begin(), lumiSections.end());
0425 
0426   edm::LogInfo("INFO") << "Number of Lumisections " << lumiSections.size() << " in run " << runSeg.run();
0427 
0428   //Determining integer number of blocks
0429   float nBlocks_f = float(lumiSections.size()) / approxLumiBlockSize_;
0430   unsigned int nBlocks = 1;
0431   if (nBlocks_f > 1) {
0432     if (nBlocks_f - lumiSections.size() / approxLumiBlockSize_ < 0.5) {
0433       nBlocks = lumiSections.size() / approxLumiBlockSize_;
0434     } else {
0435       nBlocks = lumiSections.size() / approxLumiBlockSize_ + 1;
0436     }
0437   }
0438 
0439   float nLSPerBlock = float(lumiSections.size()) / nBlocks;
0440 
0441   std::vector<std::pair<unsigned int, unsigned int>> lsKeys;
0442   lsKeys.clear();
0443 
0444   //Constructing nBlocks IOVs
0445   for (unsigned iKey = 0; iKey < nBlocks; iKey++) {
0446     lsKeys.push_back(std::make_pair(lumiSections[(unsigned int)(iKey * nLSPerBlock)],
0447                                     lumiSections[(unsigned int)((iKey + 1) * nLSPerBlock) - 1]));
0448   }
0449 
0450   lsKeys[0].first = 1;
0451 
0452   for (unsigned int lumiSection = 0; lumiSection < lumiSections.size(); lumiSection++) {
0453     thisLS = lumiSections[lumiSection];
0454 
0455     std::pair<unsigned int, unsigned int> lsKey;
0456 
0457     bool foundKey = false;
0458 
0459     for (unsigned iKey = 0; iKey < nBlocks; iKey++) {
0460       if ((thisLS >= lsKeys[iKey].first) && (thisLS <= lsKeys[iKey].second)) {
0461         lsKey = lsKeys[iKey];
0462         foundKey = true;
0463         break;
0464       }
0465     }
0466 
0467     if (!foundKey) {
0468       edm::LogInfo("WARNING") << "Didn't find key " << thisLS;
0469       continue;
0470     }
0471 
0472     if (lumiInfoMap.count(lsKey) == 0) {
0473       lumiInfoMap[lsKey] = new LumiInfo();
0474     }
0475 
0476     //Sum all lumi in IOV of lsKey
0477     totalLumiByBX_ = lumiInfoMap[lsKey]->getInstLumiAllBX();
0478     events_ = lumiInfoMap[lsKey]->getErrorLumiAllBX();
0479 
0480     rawlumiBX_ = lumiInfoMapPerLS[thisLS]->getInstLumiAllBX();
0481     errOnLumiByBX_ = lumiInfoMapPerLS[thisLS]->getErrorLumiAllBX();
0482 
0483     for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
0484       totalLumiByBX_[bx] += rawlumiBX_[bx];
0485       events_[bx] += errOnLumiByBX_[bx];
0486     }
0487     lumiInfoMap[lsKey]->setInstLumiAllBX(totalLumiByBX_);
0488     lumiInfoMap[lsKey]->setErrorLumiAllBX(events_);
0489     lumiInfoCounter[lsKey]++;
0490   }
0491 
0492   cond::Time_t thisIOV = 1;
0493 
0494   char* histname1 = new char[100];
0495   char* histname2 = new char[100];
0496   char* histname3 = new char[100];
0497   char* histTitle1 = new char[100];
0498   char* histTitle2 = new char[100];
0499   char* histTitle3 = new char[100];
0500   sprintf(histTitle1, "Type1Fraction_%d", runSeg.run());
0501   sprintf(histTitle2, "Type1Res_%d", runSeg.run());
0502   sprintf(histTitle3, "Type2Res_%d", runSeg.run());
0503   type1FracHist = new TH1F(histTitle1, histTitle1, 1000, -0.5, 0.5);
0504   type1resHist = new TH1F(histTitle2, histTitle2, 4000, -0.2, 0.2);
0505   type2resHist = new TH1F(histTitle3, histTitle3, 4000, -0.2, 0.2);
0506   delete[] histTitle1;
0507   delete[] histTitle2;
0508   delete[] histTitle3;
0509 
0510   for (lumiInfoMapIterator = lumiInfoMap.begin(); (lumiInfoMapIterator != lumiInfoMap.end()); ++lumiInfoMapIterator) {
0511     totalLumiByBX_ = lumiInfoMapIterator->second->getInstLumiAllBX();
0512     events_ = lumiInfoMapIterator->second->getErrorLumiAllBX();
0513 
0514     if (events_.empty()) {
0515       continue;
0516     }
0517 
0518     edm::LuminosityBlockID lu(runSeg.id().run(), edm::LuminosityBlockNumber_t(lumiInfoMapIterator->first.first));
0519     thisIOV = (cond::Time_t)(lu.value());
0520 
0521     sprintf(histname1,
0522             "CorrectedLumiAvg_%d_%d_%d_%d",
0523             runSeg.run(),
0524             iBlock,
0525             lumiInfoMapIterator->first.first,
0526             lumiInfoMapIterator->first.second);
0527     sprintf(histname2,
0528             "ScaleFactorsAvg_%d_%d_%d_%d",
0529             runSeg.run(),
0530             iBlock,
0531             lumiInfoMapIterator->first.first,
0532             lumiInfoMapIterator->first.second);
0533     sprintf(histname3,
0534             "RawLumiAvg_%d_%d_%d_%d",
0535             runSeg.run(),
0536             iBlock,
0537             lumiInfoMapIterator->first.first,
0538             lumiInfoMapIterator->first.second);
0539 
0540     corrlumiAvg_h = new TH1F(histname1, "", LumiConstants::numBX, 1, LumiConstants::numBX);
0541     scaleFactorAvg_h = new TH1F(histname2, "", LumiConstants::numBX, 1, LumiConstants::numBX);
0542     lumiAvg_h = new TH1F(histname3, "", LumiConstants::numBX, 1, LumiConstants::numBX);
0543 
0544     //Averaging by the number of events
0545     for (unsigned int i = 0; i < LumiConstants::numBX; i++) {
0546       if (events_.at(i) != 0) {
0547         totalLumiByBXAvg_[i] = totalLumiByBX_[i] / events_[i];
0548       } else {
0549         totalLumiByBXAvg_[i] = 0.0;
0550       }
0551     }
0552 
0553     calculateCorrections(totalLumiByBXAvg_, correctionScaleFactors_, overallCorrection_);
0554 
0555     for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
0556       corrlumiAvg_h->SetBinContent(bx, totalLumiByBXAvg_[bx] * correctionScaleFactors_[bx]);
0557       if (events_.at(bx) != 0) {
0558         corrlumiAvg_h->SetBinError(bx,
0559                                    totalLumiByBXAvg_[bx] * correctionScaleFactors_[bx] / TMath::Sqrt(events_.at(bx)));
0560       } else {
0561         corrlumiAvg_h->SetBinError(bx, 0.0);
0562       }
0563 
0564       scaleFactorAvg_h->SetBinContent(bx, correctionScaleFactors_[bx]);
0565       lumiAvg_h->SetBinContent(bx, totalLumiByBXAvg_[bx]);
0566     }
0567 
0568     //Writing the corrections to SQL lite file for db
0569     LumiCorrections pccCorrections;
0570     pccCorrections.setOverallCorrection(overallCorrection_);
0571     pccCorrections.setType1Fraction(type1Frac);
0572     pccCorrections.setType1Residual(mean_type1_residual);
0573     pccCorrections.setType2Residual(mean_type2_residual);
0574     pccCorrections.setCorrectionsBX(correctionScaleFactors_);
0575 
0576     if (poolDbService.isAvailable()) {
0577       poolDbService->writeOneIOV(pccCorrections, thisIOV, "LumiCorrectionsRcd");
0578     } else {
0579       throw std::runtime_error("PoolDBService required.");
0580     }
0581 
0582     histoFile->cd();
0583     corrlumiAvg_h->Write();
0584     scaleFactorAvg_h->Write();
0585     lumiAvg_h->Write();
0586 
0587     delete corrlumiAvg_h;
0588     delete scaleFactorAvg_h;
0589     delete lumiAvg_h;
0590 
0591     type1FracHist->Fill(type1Frac);
0592     type1resHist->Fill(mean_type1_residual);
0593     type2resHist->Fill(mean_type2_residual);
0594 
0595     for (unsigned int ils = lumiInfoMapIterator->first.first; ils < lumiInfoMapIterator->first.second + 1; ils++) {
0596       if (ils > maxLS) {
0597         std::cout << "ils out of maxLS range!!" << std::endl;
0598         break;
0599       }
0600       Type1FracMon->setBinContent(ils, type1Frac);
0601       Type1FracMon->setBinError(ils, mean_type1_residual_unc);
0602       Type1ResMon->setBinContent(ils, mean_type1_residual);
0603       Type1ResMon->setBinError(ils, mean_type1_residual_unc);
0604       Type2ResMon->setBinContent(ils, mean_type2_residual);
0605       Type2ResMon->setBinError(ils, mean_type2_residual_unc);
0606     }
0607 
0608     type1FracGraph->SetPoint(iBlock, thisIOV + approxLumiBlockSize_ / 2.0, type1Frac);
0609     type1resGraph->SetPoint(iBlock, thisIOV + approxLumiBlockSize_ / 2.0, mean_type1_residual);
0610     type2resGraph->SetPoint(iBlock, thisIOV + approxLumiBlockSize_ / 2.0, mean_type2_residual);
0611     type1FracGraph->SetPointError(iBlock, approxLumiBlockSize_ / 2.0, mean_type1_residual_unc);
0612     type1resGraph->SetPointError(iBlock, approxLumiBlockSize_ / 2.0, mean_type1_residual_unc);
0613     type2resGraph->SetPointError(iBlock, approxLumiBlockSize_ / 2.0, mean_type2_residual_unc);
0614     pedestalGraph->SetPoint(iBlock, thisIOV + approxLumiBlockSize_ / 2.0, pedestal);
0615     pedestalGraph->SetPointError(iBlock, approxLumiBlockSize_ / 2.0, pedestal_unc);
0616 
0617     edm::LogInfo("INFO")
0618         << "iBlock type1Frac mean_type1_residual mean_type2_residual mean_type1_residual_unc mean_type2_residual_unc "
0619         << iBlock << " " << type1Frac << " " << mean_type1_residual << " " << mean_type2_residual << " "
0620         << mean_type1_residual_unc << " " << mean_type2_residual_unc;
0621 
0622     type1Frac = 0.0;
0623     mean_type1_residual = 0.0;
0624     mean_type2_residual = 0.0;
0625     mean_type1_residual_unc = 0;
0626     mean_type2_residual_unc = 0;
0627     pedestal = 0.;
0628     pedestal_unc = 0.;
0629 
0630     iBlock++;
0631 
0632     resetBlock();
0633   }
0634   histoFile->cd();
0635   type1FracHist->Write();
0636   type1resHist->Write();
0637   type2resHist->Write();
0638 
0639   delete type1FracHist;
0640   delete type1resHist;
0641   delete type2resHist;
0642 
0643   delete[] histname1;
0644   delete[] histname2;
0645   delete[] histname3;
0646 
0647   for (lumiInfoMapIterator = lumiInfoMap.begin(); (lumiInfoMapIterator != lumiInfoMap.end()); ++lumiInfoMapIterator) {
0648     delete lumiInfoMapIterator->second;
0649   }
0650   for (unsigned int lumiSection = 0; lumiSection < lumiSections.size(); lumiSection++) {
0651     thisLS = lumiSections[lumiSection];
0652     delete lumiInfoMapPerLS[thisLS];
0653   }
0654   lumiInfoMap.clear();
0655   lumiInfoMapPerLS.clear();
0656   lumiSections.clear();
0657   lumiInfoCounter.clear();
0658 }
0659 
0660 //--------------------------------------------------------------------------------------------------
0661 void CorrPCCProducer::resetBlock() {
0662   for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
0663     totalLumiByBX_[bx] = 0;
0664     totalLumiByBXAvg_[bx] = 0;
0665     events_[bx] = 0;
0666     correctionScaleFactors_[bx] = 1.0;
0667   }
0668 }
0669 //--------------------------------------------------------------------------------------------------
0670 void CorrPCCProducer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& context) {
0671   ibooker.setCurrentFolder("AlCaReco/LumiPCC/");
0672   auto scope = DQMStore::IBooker::UseRunScope(ibooker);
0673   Type1FracMon = ibooker.book1D("type1Fraction", "Type1Fraction;Lumisection;Type 1 Fraction", maxLS, 0, maxLS);
0674   Type1ResMon = ibooker.book1D("type1Residual", "Type1Residual;Lumisection;Type 1 Residual", maxLS, 0, maxLS);
0675   Type2ResMon = ibooker.book1D("type2Residual", "Type2Residual;Lumisection;Type 2 Residual", maxLS, 0, maxLS);
0676 }
0677 //--------------------------------------------------------------------------------------------------
0678 void CorrPCCProducer::endJob() {
0679   histoFile->cd();
0680   type1FracGraph->Write();
0681   type1resGraph->Write();
0682   type2resGraph->Write();
0683   pedestalGraph->Write();
0684   histoFile->Write();
0685   histoFile->Close();
0686   delete type1FracGraph;
0687   delete type1resGraph;
0688   delete type2resGraph;
0689   delete pedestalGraph;
0690 }
0691 
0692 DEFINE_FWK_MODULE(CorrPCCProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!