14 #include <boost/serialization/vector.hpp> 45 #include "TGraphErrors.h" 95 std::map<std::pair<unsigned int, unsigned int>,
LumiInfo*>
160 lumiInfoToken = consumes<LumiInfo, edm::InLumi>(inputPCCTag_);
162 histoFile =
new TFile(
"CorrectionHisto.root",
"RECREATE");
196 for (
size_t i = 0;
i < lumi_vector.size();
i++) {
197 if (lumi_vector.at(
i) > max_lumi)
198 max_lumi = lumi_vector.at(
i);
207 std::vector<float> corrected_tmp_;
208 for (
size_t i = 0;
i < uncorrPCCPerBX.size();
i++) {
209 corrected_tmp_.push_back(uncorrPCCPerBX.at(
i));
214 float bin_k = corrected_tmp_.at(
k);
215 corrected_tmp_.at(
k + 1) = corrected_tmp_.at(
k + 1) - type1Frac * bin_k;
219 for (
size_t i = 0;
i < LumiConstants::numBX - 1;
i++) {
220 for (
size_t j =
i + 1; j <
i + LumiConstants::numBX - 1; j++) {
221 float bin_i = corrected_tmp_.at(
i);
222 if (j < LumiConstants::numBX) {
225 corrected_tmp_.at(j - LumiConstants::numBX) =
245 std::vector<float> afterGlow;
246 TH1F type1(
"type1",
"", 1000, -0.5, 0.5);
247 TH1F type2(
"type2",
"", 1000, -0.5, 0.5);
249 lumi = corrected_tmp_.at(ibx);
251 afterGlow.push_back(corrected_tmp_.at(ibx + 1));
252 afterGlow.push_back(corrected_tmp_.at(ibx + 2));
255 if (lumi > threshold && afterGlow[0] < threshold && afterGlow[1] < threshold) {
257 float thisAfterGlow = corrected_tmp_.at(ibx +
index);
258 if (thisAfterGlow < threshold) {
259 afterGlow.push_back(thisAfterGlow);
266 if (afterGlow.size() > 1) {
269 thisType2 += afterGlow[
index];
270 type2.Fill(afterGlow[
index] / lumi);
275 thisType1 = (afterGlow[0] - thisType2) / lumi;
277 type1.Fill(thisType1);
300 for (
int trial = 0; trial < nTrials; trial++) {
308 std::vector<float> corrected_tmp_;
309 for (
size_t i = 0;
i < uncorrected.size();
i++) {
310 corrected_tmp_.push_back(uncorrected.at(
i));
316 float bin_i = corrected_tmp_.at(
i);
317 corrected_tmp_.at(
i + 1) = corrected_tmp_.at(
i + 1) - type1Frac * bin_i;
320 bin_i = corrected_tmp_.at(
i);
321 for (
size_t j =
i + 1; j <
i + LumiConstants::numBX - 1; j++) {
322 if (j < LumiConstants::numBX) {
325 corrected_tmp_.at(j - LumiConstants::numBX) =
333 float integral_uncorr_clusters = 0;
334 float integral_corr_clusters = 0;
339 for (
size_t ibx = 0; ibx < corrected_tmp_.size(); ibx++) {
340 if (corrected_tmp_.at(ibx) >
threshold) {
341 integral_uncorr_clusters += uncorrected.at(ibx);
342 integral_corr_clusters += corrected_tmp_.at(ibx);
344 if (corrected_tmp_.at(ibx) != 0.0 && uncorrected.at(ibx) != 0.0) {
345 correctionScaleFactors_.at(ibx) = corrected_tmp_.at(ibx) / uncorrected.at(ibx);
347 correctionScaleFactors_.at(ibx) = 0.0;
351 overallCorrection_ = integral_corr_clusters / integral_uncorr_clusters;
370 unsigned int totalEvents = 0;
376 edm::LogInfo(
"INFO") <<
"number of events in this LS is too few " << totalEvents;
405 unsigned int nBlocks = 1;
416 std::vector<std::pair<unsigned int, unsigned int>> lsKeys;
420 for (
unsigned iKey = 0; iKey < nBlocks; iKey++) {
421 lsKeys.push_back(std::make_pair(
lumiSections[(
unsigned int)(iKey * nLSPerBlock)],
422 lumiSections[(
unsigned int)((iKey + 1) * nLSPerBlock) - 1]));
427 for (
unsigned int lumiSection = 0; lumiSection <
lumiSections.size(); lumiSection++) {
430 std::pair<unsigned int, unsigned int> lsKey;
432 bool foundKey =
false;
434 for (
unsigned iKey = 0; iKey < nBlocks; iKey++) {
436 lsKey = lsKeys[iKey];
469 char* histname1 =
new char[100];
470 char* histname2 =
new char[100];
471 char* histname3 =
new char[100];
472 char* histTitle1 =
new char[100];
473 char* histTitle2 =
new char[100];
474 char* histTitle3 =
new char[100];
475 sprintf(histTitle1,
"Type1Fraction_%d", runSeg.
run());
476 sprintf(histTitle2,
"Type1Res_%d", runSeg.
run());
477 sprintf(histTitle3,
"Type2Res_%d", runSeg.
run());
478 type1FracHist =
new TH1F(histTitle1, histTitle1, 1000, -0.5, 0.5);
479 type1resHist =
new TH1F(histTitle2, histTitle2, 4000, -0.2, 0.2);
480 type2resHist =
new TH1F(histTitle3, histTitle3, 4000, -0.2, 0.2);
497 "CorrectedLumiAvg_%d_%d_%d_%d",
503 "ScaleFactorsAvg_%d_%d_%d_%d",
509 "RawLumiAvg_%d_%d_%d_%d",
533 corrlumiAvg_h->SetBinError(bx,
536 corrlumiAvg_h->SetBinError(bx, 0.0);
554 throw std::runtime_error(
"PoolDBService required.");
560 corrlumiAvg_h->Write();
574 std::cout <<
"ils out of maxLS range!!" << std::endl;
593 <<
"iBlock type1Frac mean_type1_residual mean_type2_residual mean_type1_residual_unc mean_type2_residual_unc " 601 mean_type2_residual_unc = 0;
608 type1FracHist->Write();
623 for (
unsigned int lumiSection = 0; lumiSection <
lumiSections.size(); lumiSection++) {
T getParameter(std::string const &) const
std::vector< float > events_
void calculateCorrections(std::vector< float >, std::vector< float > &, float &)
MonitorElement * Type1FracMon
void setBinContent(int binx, double content)
set content of bin (1-D)
std::map< std::pair< unsigned int, unsigned int >, LumiInfo * > lumiInfoMap
void setType1Fraction(float type1frac)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
float mean_type1_residual
edm::EDGetTokenT< LumiInfo > lumiInfoToken
void endRun(edm::Run const &runSeg, const edm::EventSetup &iSetup) final
void setOverallCorrection(float overallCorrection)
static const unsigned int numBX
const std::vector< float > & getErrorLumiAllBX() const
unsigned int LuminosityBlockNumber_t
MonitorElement * Type1ResMon
LuminosityBlockNumber_t luminosityBlock() const
void beginLuminosityBlock(edm::LuminosityBlock const &lumiSeg, const edm::EventSetup &iSetup) final
#define DEFINE_FWK_MODULE(type)
unsigned long long Time_t
void endRunProduce(edm::Run &runSeg, const edm::EventSetup &iSetup) final
void setCurrentFolder(std::string const &fullpath)
float mean_type1_residual_unc
unsigned int minimumNumberOfEvents
std::vector< float > totalLumiByBXAvg_
std::map< std::pair< unsigned int, unsigned int >, unsigned int > lumiInfoCounter
void evaluateCorrectionResiduals(std::vector< float >)
TGraphErrors * type1resGraph
LumiCorrections * pccCorrections
MonitorElement * book1D(Args &&...args)
void writeOne(T *payload, Time_t time, const std::string &recordName, bool withlogging=false)
void setBinError(int binx, double error)
set uncertainty on content of bin (1-D)
std::vector< float > correctionScaleFactors_
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
void setType2Residual(float type2res)
const std::vector< float > & getInstLumiAllBX() const
std::vector< float > correctionTemplate_
std::vector< unsigned int > lumiSections
float getMaximum(std::vector< float >)
TGraphErrors * type2resGraph
std::vector< float > rawlumiBX_
T const * product() const
std::vector< float > totalLumiByBX_
float mean_type2_residual_unc
edm::Service< cond::service::PoolDBOutputService > poolDbService
TGraphErrors * type1FracGraph
unsigned int approxLumiBlockSize_
void endLuminosityBlock(edm::LuminosityBlock const &lumiSeg, const edm::EventSetup &iSetup) final
void makeCorrectionTemplate()
void estimateType1Frac(std::vector< float >, float &)
~CorrPCCProducer() override
std::map< std::pair< unsigned int, unsigned int >, LumiInfo * >::iterator lumiInfoMapIterator
void setType1Residual(float type1res)
float mean_type2_residual
CorrPCCProducer(const edm::ParameterSet &)
std::vector< float > errOnLumiByBX_
MonitorElement * Type2ResMon
void setCorrectionsBX(std::vector< float > &correctBX)
std::map< unsigned int, LumiInfo * > lumiInfoMapPerLS