Inheritance diagram for CorrPCCProducer:

Public Member Functions
	CorrPCCProducer (const edm::ParameterSet &)

	~CorrPCCProducer () override

Public Member Functions inherited from one::DQMEDAnalyzer< edm::one::WatchLuminosityBlocks >
	DQMEDAnalyzer ()=default

	DQMEDAnalyzer (DQMEDAnalyzer< T... > const &)=delete

	DQMEDAnalyzer (DQMEDAnalyzer< T... > &&)=delete

	~DQMEDAnalyzer () override=default

Private Member Functions
void	beginLuminosityBlock (edm::LuminosityBlock const &lumiSeg, const edm::EventSetup &iSetup) final

void	bookHistograms (DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override

void	calculateCorrections (std::vector< float >, std::vector< float > &, float &)

void	endJob () final

void	endLuminosityBlock (edm::LuminosityBlock const &lumiSeg, const edm::EventSetup &iSetup) final

void	endRun (edm::Run const &runSeg, const edm::EventSetup &iSetup) final

void	endRunProduce (edm::Run &runSeg, const edm::EventSetup &iSetup) final

void	estimateType1Frac (std::vector< float >, float &)

void	evaluateCorrectionResiduals (std::vector< float >)

float	getMaximum (std::vector< float >)

void	makeCorrectionTemplate ()

void	resetBlock ()

Private Attributes
unsigned int	approxLumiBlockSize_

std::vector< float >	correctionScaleFactors_

std::vector< float >	correctionTemplate_

TH1F *	corrlumiAvg_h

unsigned int	countLumi_

std::vector< float >	errOnLumiByBX_

std::vector< float >	events_

TFile *	histoFile

TList *	hlist

unsigned int	iBlock = 0

TH1F *	lumiAvg_h

std::map< std::pair< unsigned int, unsigned int >, unsigned int >	lumiInfoCounter

std::map< std::pair< unsigned int, unsigned int >, LumiInfo * >	lumiInfoMap

std::map< std::pair< unsigned int, unsigned int >, LumiInfo * >::iterator	lumiInfoMapIterator

std::map< unsigned int, LumiInfo * >	lumiInfoMapPerLS

edm::EDGetTokenT< LumiInfo >	lumiInfoToken

std::vector< unsigned int >	lumiSections

unsigned int	maxLS = 3500

float	mean_type1_residual

float	mean_type1_residual_unc

float	mean_type2_residual

float	mean_type2_residual_unc

unsigned int	minimumNumberOfEvents

unsigned int	nTrain

float	overallCorrection_

LumiCorrections *	pccCorrections

std::string	pccSrc_

edm::Service< cond::service::PoolDBOutputService >	poolDbService

std::string	prodInst_

std::vector< float >	rawlumiBX_

TH1F *	scaleFactorAvg_h

unsigned int	thisLS

std::vector< float >	totalLumiByBX_

std::vector< float >	totalLumiByBXAvg_

float	type1Frac

TGraphErrors *	type1FracGraph

TH1F *	type1FracHist

MonitorElement *	Type1FracMon

TGraphErrors *	type1resGraph

TH1F *	type1resHist

MonitorElement *	Type1ResMon

double	type2_a_

double	type2_b_

TGraphErrors *	type2resGraph

TH1F *	type2resHist

MonitorElement *	Type2ResMon

Detailed Description

Definition at line 55 of file CorrPCCProducer.cc.

Constructor & Destructor Documentation

CorrPCCProducer::CorrPCCProducer ( const edm::ParameterSet & iConfig )

explicit

Definition at line 136 of file CorrPCCProducer.cc.

References approxLumiBlockSize_, correctionScaleFactors_, correctionTemplate_, countLumi_, events_, edm::ParameterSet::getParameter(), histoFile, lumiInfoToken, makeCorrectionTemplate(), minimumNumberOfEvents, LumiConstants::numBX, pccSrc_, prodInst_, resetBlock(), totalLumiByBX_, totalLumiByBXAvg_, type1FracGraph, type1resGraph, type2_a_, type2_b_, and type2resGraph.

                                                                {
   pccSrc_ =
       iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<std::string>("inLumiObLabel");
   prodInst_ =
       iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<std::string>("ProdInst");
   approxLumiBlockSize_ =
       iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<int>("approxLumiBlockSize");
   type2_a_ = iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<double>("type2_a");
   type2_b_ = iConfig.getParameter<edm::ParameterSet>("CorrPCCProducerParameters").getParameter<double>("type2_b");
   countLumi_ = 0;
   minimumNumberOfEvents = 1000;
 
   totalLumiByBX_.resize(LumiConstants::numBX);
   totalLumiByBXAvg_.resize(LumiConstants::numBX);
   events_.resize(LumiConstants::numBX);
   correctionScaleFactors_.resize(LumiConstants::numBX);
   correctionTemplate_.resize(LumiConstants::numBX);
 
   resetBlock();
 
   makeCorrectionTemplate();
 
   edm::InputTag inputPCCTag_(pccSrc_, prodInst_);
 
   lumiInfoToken = consumes<LumiInfo, edm::InLumi>(inputPCCTag_);
 
   histoFile = new TFile("CorrectionHisto.root", "RECREATE");
 
   type1FracGraph = new TGraphErrors();
   type1resGraph = new TGraphErrors();
   type2resGraph = new TGraphErrors();
   type1FracGraph->SetName("Type1Fraction");
   type1resGraph->SetName("Type1Res");
   type2resGraph->SetName("Type2Res");
   type1FracGraph->GetYaxis()->SetTitle("Type 1 Fraction");
   type1resGraph->GetYaxis()->SetTitle("Type 1 Residual");
   type2resGraph->GetYaxis()->SetTitle("Type 2 Residual");
   type1FracGraph->GetXaxis()->SetTitle("Unique LS ID");
   type1resGraph->GetXaxis()->SetTitle("Unique LS ID");
   type2resGraph->GetXaxis()->SetTitle("Unique LS ID");
   type1FracGraph->SetMarkerStyle(8);
   type1resGraph->SetMarkerStyle(8);
   type2resGraph->SetMarkerStyle(8);
 }

CorrPCCProducer::~CorrPCCProducer ( )

override

Definition at line 182 of file CorrPCCProducer.cc.

182 {}

Member Function Documentation

void CorrPCCProducer::beginLuminosityBlock	(	edm::LuminosityBlock const &	lumiSeg,
		const edm::EventSetup &	iSetup
	)

finalprivate

Definition at line 355 of file CorrPCCProducer.cc.

References countLumi_.

                                                                                                          {
   countLumi_++;
 }

void CorrPCCProducer::bookHistograms	(	DQMStore::IBooker &	ibooker,
		edm::Run const &	iRun,
		edm::EventSetup const &	context
	)

overrideprivate

Definition at line 643 of file CorrPCCProducer.cc.

References DQMStore::IBooker::book1D(), maxLS, DQMStore::IBooker::setCurrentFolder(), Type1FracMon, Type1ResMon, and Type2ResMon.

                                                                                                                {
   ibooker.setCurrentFolder("AlCaReco/LumiPCC/");
   Type1FracMon = ibooker.book1D("type1Fraction", "Type1Fraction;Lumisection;Type 1 Fraction", maxLS, 0, maxLS);
   Type1ResMon = ibooker.book1D("type1Residual", "Type1Residual;Lumisection;Type 1 Residual", maxLS, 0, maxLS);
   Type2ResMon = ibooker.book1D("type2Residual", "Type2Residual;Lumisection;Type 2 Residual", maxLS, 0, maxLS);
 }

void CorrPCCProducer::calculateCorrections	(	std::vector< float >	uncorrected,
		std::vector< float > &	correctionScaleFactors_,
		float &	overallCorrection_
	)

private

Definition at line 293 of file CorrPCCProducer.cc.

References correctionTemplate_, estimateType1Frac(), evaluateCorrectionResiduals(), getMaximum(), mps_fire::i, SiStripPI::max, LumiConstants::numBX, electronIdCutBased_cfi::threshold, and type1Frac.

Referenced by endRunProduce().

                                                                       {
   type1Frac = 0;
 
   int nTrials = 4;
 
   for (int trial = 0; trial < nTrials; trial++) {
     estimateType1Frac(uncorrected, type1Frac);
     edm::LogInfo("INFO") << "type 1 fraction after iteration " << trial << " is  " << type1Frac;
   }
 
   //correction should never be negative
   type1Frac = std::max(0.0, (double)type1Frac);
 
   std::vector<float> corrected_tmp_;
   for (size_t i = 0; i < uncorrected.size(); i++) {
     corrected_tmp_.push_back(uncorrected.at(i));
   }
 
   //Apply all corrections
   for (size_t i = 0; i < LumiConstants::numBX - 1; i++) {
     // type 1 - first (spillover from previous BXs real clusters)
     float bin_i = corrected_tmp_.at(i);
     corrected_tmp_.at(i + 1) = corrected_tmp_.at(i + 1) - type1Frac * bin_i;
 
     // type 2 - after (comes from real activation)
     bin_i = corrected_tmp_.at(i);
     for (size_t j = i + 1; j < i + LumiConstants::numBX - 1; j++) {
       if (j < LumiConstants::numBX) {
         corrected_tmp_.at(j) = corrected_tmp_.at(j) - bin_i * correctionTemplate_.at(j - i);
       } else {
         corrected_tmp_.at(j - LumiConstants::numBX) =
             corrected_tmp_.at(j - LumiConstants::numBX) - bin_i * correctionTemplate_.at(j - i);
       }
     }
   }
 
   evaluateCorrectionResiduals(corrected_tmp_);
 
   float integral_uncorr_clusters = 0;
   float integral_corr_clusters = 0;
 
   //Calculate Per-BX correction factor and overall correction factor
   float lumiMax = getMaximum(corrected_tmp_);
   float threshold = lumiMax * 0.2;
   for (size_t ibx = 0; ibx < corrected_tmp_.size(); ibx++) {
     if (corrected_tmp_.at(ibx) > threshold) {
       integral_uncorr_clusters += uncorrected.at(ibx);
       integral_corr_clusters += corrected_tmp_.at(ibx);
     }
     if (corrected_tmp_.at(ibx) != 0.0 && uncorrected.at(ibx) != 0.0) {
       correctionScaleFactors_.at(ibx) = corrected_tmp_.at(ibx) / uncorrected.at(ibx);
     } else {
       correctionScaleFactors_.at(ibx) = 0.0;
     }
   }
 
   overallCorrection_ = integral_corr_clusters / integral_uncorr_clusters;
 }

void CorrPCCProducer::endJob ( void )

finalprivate

Definition at line 650 of file CorrPCCProducer.cc.

References DEFINE_FWK_MODULE, histoFile, type1FracGraph, type1resGraph, and type2resGraph.

Referenced by o2olib.O2ORunMgr::executeJob().

                              {
   histoFile->cd();
   type1FracGraph->Write();
   type1resGraph->Write();
   type2resGraph->Write();
   histoFile->Write();
   histoFile->Close();
   delete type1FracGraph;
   delete type1resGraph;
   delete type2resGraph;
 }

void CorrPCCProducer::endLuminosityBlock	(	edm::LuminosityBlock const &	lumiSeg,
		const edm::EventSetup &	iSetup
	)

finalprivate

Definition at line 360 of file CorrPCCProducer.cc.

References errOnLumiByBX_, events_, edm::LuminosityBlock::getByToken(), LumiInfo::getErrorLumiAllBX(), LumiInfo::getInstLumiAllBX(), lumiInfoMapPerLS, lumiInfoToken, edm::LuminosityBlockBase::luminosityBlock(), lumiSections, minimumNumberOfEvents, LumiConstants::numBX, edm::Handle< T >::product(), thisLS, and totalLumiByBX_.

                                                                                                        {
   thisLS = lumiSeg.luminosityBlock();
 
   edm::Handle<LumiInfo> PCCHandle;
   lumiSeg.getByToken(lumiInfoToken, PCCHandle);
 
   const LumiInfo& inLumiOb = *(PCCHandle.product());
 
   errOnLumiByBX_ = inLumiOb.getErrorLumiAllBX();
 
   unsigned int totalEvents = 0;
   for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
     totalEvents += errOnLumiByBX_[bx];
   }
 
   if (totalEvents < minimumNumberOfEvents) {
     edm::LogInfo("INFO") << "number of events in this LS is too few " << totalEvents;
     //return;
   } else {
     edm::LogInfo("INFO") << "Skipping Lumisection " << thisLS;
   }
 
   lumiInfoMapPerLS[thisLS] = new LumiInfo();
   totalLumiByBX_ = inLumiOb.getInstLumiAllBX();
   events_ = inLumiOb.getErrorLumiAllBX();
   lumiInfoMapPerLS[thisLS]->setInstLumiAllBX(totalLumiByBX_);
   lumiInfoMapPerLS[thisLS]->setErrorLumiAllBX(events_);
   lumiSections.push_back(thisLS);
 }

void CorrPCCProducer::endRun	(	edm::Run const &	runSeg,
		const edm::EventSetup &	iSetup
	)

finalprivate

Definition at line 391 of file CorrPCCProducer.cc.

391 {}

void CorrPCCProducer::endRunProduce	(	edm::Run &	runSeg,
		const edm::EventSetup &	iSetup
	)

finalprivate

Definition at line 394 of file CorrPCCProducer.cc.

References approxLumiBlockSize_, calculateCorrections(), correctionScaleFactors_, corrlumiAvg_h, gather_cfg::cout, errOnLumiByBX_, events_, plotBeamSpotDB::first, objects.autophobj::float, histoFile, mps_fire::i, iBlock, edm::RunBase::id(), edm::Service< T >::isAvailable(), lumiAvg_h, lumiInfoCounter, lumiInfoMap, lumiInfoMapIterator, lumiInfoMapPerLS, lumiSections, maxLS, mean_type1_residual, mean_type1_residual_unc, mean_type2_residual, mean_type2_residual_unc, LumiConstants::numBX, overallCorrection_, pccCorrections, poolDbService, rawlumiBX_, resetBlock(), edm::RunID::run(), edm::RunBase::run(), scaleFactorAvg_h, MonitorElement::setBinContent(), MonitorElement::setBinError(), LumiCorrections::setCorrectionsBX(), LumiCorrections::setOverallCorrection(), LumiCorrections::setType1Fraction(), LumiCorrections::setType1Residual(), LumiCorrections::setType2Residual(), jetUpdater_cfi::sort, thisLS, totalLumiByBX_, totalLumiByBXAvg_, type1Frac, type1FracGraph, type1FracHist, Type1FracMon, type1resGraph, type1resHist, Type1ResMon, type2resGraph, type2resHist, Type2ResMon, and cond::service::PoolDBOutputService::writeOne().

                                                                                {
   if (lumiSections.empty()) {
     return;
   }
 
   std::sort(lumiSections.begin(), lumiSections.end());
 
   edm::LogInfo("INFO") << "Number of Lumisections " << lumiSections.size() << " in run " << runSeg.run();
 
   //Determining integer number of blocks
   float nBlocks_f = float(lumiSections.size()) / approxLumiBlockSize_;
   unsigned int nBlocks = 1;
   if (nBlocks_f > 1) {
     if (nBlocks_f - lumiSections.size() / approxLumiBlockSize_ < 0.5) {
       nBlocks = lumiSections.size() / approxLumiBlockSize_;
     } else {
       nBlocks = lumiSections.size() / approxLumiBlockSize_ + 1;
     }
   }
 
   float nLSPerBlock = float(lumiSections.size()) / nBlocks;
 
   std::vector<std::pair<unsigned int, unsigned int>> lsKeys;
   lsKeys.clear();
 
   //Constructing nBlocks IOVs
   for (unsigned iKey = 0; iKey < nBlocks; iKey++) {
     lsKeys.push_back(std::make_pair(lumiSections[(unsigned int)(iKey * nLSPerBlock)],
                                     lumiSections[(unsigned int)((iKey + 1) * nLSPerBlock) - 1]));
   }
 
   lsKeys[0].first = 1;
 
   for (unsigned int lumiSection = 0; lumiSection < lumiSections.size(); lumiSection++) {
     thisLS = lumiSections[lumiSection];
 
     std::pair<unsigned int, unsigned int> lsKey;
 
     bool foundKey = false;
 
     for (unsigned iKey = 0; iKey < nBlocks; iKey++) {
       if ((thisLS >= lsKeys[iKey].first) && (thisLS <= lsKeys[iKey].second)) {
         lsKey = lsKeys[iKey];
         foundKey = true;
         break;
       }
     }
 
     if (!foundKey) {
       edm::LogInfo("WARNING") << "Didn't find key " << thisLS;
       continue;
     }
 
     if (lumiInfoMap.count(lsKey) == 0) {
       lumiInfoMap[lsKey] = new LumiInfo();
     }
 
     //Sum all lumi in IOV of lsKey
     totalLumiByBX_ = lumiInfoMap[lsKey]->getInstLumiAllBX();
     events_ = lumiInfoMap[lsKey]->getErrorLumiAllBX();
 
     rawlumiBX_ = lumiInfoMapPerLS[thisLS]->getInstLumiAllBX();
     errOnLumiByBX_ = lumiInfoMapPerLS[thisLS]->getErrorLumiAllBX();
 
     for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
       totalLumiByBX_[bx] += rawlumiBX_[bx];
       events_[bx] += errOnLumiByBX_[bx];
     }
     lumiInfoMap[lsKey]->setInstLumiAllBX(totalLumiByBX_);
     lumiInfoMap[lsKey]->setErrorLumiAllBX(events_);
     lumiInfoCounter[lsKey]++;
   }
 
   cond::Time_t thisIOV = 1;
 
   char* histname1 = new char[100];
   char* histname2 = new char[100];
   char* histname3 = new char[100];
   char* histTitle1 = new char[100];
   char* histTitle2 = new char[100];
   char* histTitle3 = new char[100];
   sprintf(histTitle1, "Type1Fraction_%d", runSeg.run());
   sprintf(histTitle2, "Type1Res_%d", runSeg.run());
   sprintf(histTitle3, "Type2Res_%d", runSeg.run());
   type1FracHist = new TH1F(histTitle1, histTitle1, 1000, -0.5, 0.5);
   type1resHist = new TH1F(histTitle2, histTitle2, 4000, -0.2, 0.2);
   type2resHist = new TH1F(histTitle3, histTitle3, 4000, -0.2, 0.2);
   delete[] histTitle1;
   delete[] histTitle2;
   delete[] histTitle3;
 
   for (lumiInfoMapIterator = lumiInfoMap.begin(); (lumiInfoMapIterator != lumiInfoMap.end()); ++lumiInfoMapIterator) {
     totalLumiByBX_ = lumiInfoMapIterator->second->getInstLumiAllBX();
     events_ = lumiInfoMapIterator->second->getErrorLumiAllBX();
 
     if (events_.empty()) {
       continue;
     }
 
     edm::LuminosityBlockID lu(runSeg.id().run(), edm::LuminosityBlockNumber_t(lumiInfoMapIterator->first.first));
     thisIOV = (cond::Time_t)(lu.value());
 
     sprintf(histname1,
             "CorrectedLumiAvg_%d_%d_%d_%d",
             runSeg.run(),
             iBlock,
             lumiInfoMapIterator->first.first,
             lumiInfoMapIterator->first.second);
     sprintf(histname2,
             "ScaleFactorsAvg_%d_%d_%d_%d",
             runSeg.run(),
             iBlock,
             lumiInfoMapIterator->first.first,
             lumiInfoMapIterator->first.second);
     sprintf(histname3,
             "RawLumiAvg_%d_%d_%d_%d",
             runSeg.run(),
             iBlock,
             lumiInfoMapIterator->first.first,
             lumiInfoMapIterator->first.second);
 
     corrlumiAvg_h = new TH1F(histname1, "", LumiConstants::numBX, 1, LumiConstants::numBX);
     scaleFactorAvg_h = new TH1F(histname2, "", LumiConstants::numBX, 1, LumiConstants::numBX);
     lumiAvg_h = new TH1F(histname3, "", LumiConstants::numBX, 1, LumiConstants::numBX);
 
     //Averaging by the number of events
     for (unsigned int i = 0; i < LumiConstants::numBX; i++) {
       if (events_.at(i) != 0) {
         totalLumiByBXAvg_[i] = totalLumiByBX_[i] / events_[i];
       } else {
         totalLumiByBXAvg_[i] = 0.0;
       }
     }
 
     calculateCorrections(totalLumiByBXAvg_, correctionScaleFactors_, overallCorrection_);
 
     for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
       corrlumiAvg_h->SetBinContent(bx, totalLumiByBXAvg_[bx] * correctionScaleFactors_[bx]);
       if (events_.at(bx) != 0) {
         corrlumiAvg_h->SetBinError(bx,
                                    totalLumiByBXAvg_[bx] * correctionScaleFactors_[bx] / TMath::Sqrt(events_.at(bx)));
       } else {
         corrlumiAvg_h->SetBinError(bx, 0.0);
       }
 
       scaleFactorAvg_h->SetBinContent(bx, correctionScaleFactors_[bx]);
       lumiAvg_h->SetBinContent(bx, totalLumiByBXAvg_[bx]);
     }
 
     //Writing the corrections to SQL lite file for db
     pccCorrections = new LumiCorrections();
     pccCorrections->setOverallCorrection(overallCorrection_);
     pccCorrections->setType1Fraction(type1Frac);
     pccCorrections->setType1Residual(mean_type1_residual);
     pccCorrections->setType2Residual(mean_type2_residual);
     pccCorrections->setCorrectionsBX(correctionScaleFactors_);
 
     if (poolDbService.isAvailable()) {
       poolDbService->writeOne<LumiCorrections>(pccCorrections, thisIOV, "LumiCorrectionsRcd");
     } else {
       throw std::runtime_error("PoolDBService required.");
     }
 
     delete pccCorrections;
 
     histoFile->cd();
     corrlumiAvg_h->Write();
     scaleFactorAvg_h->Write();
     lumiAvg_h->Write();
 
     delete corrlumiAvg_h;
     delete scaleFactorAvg_h;
     delete lumiAvg_h;
 
     type1FracHist->Fill(type1Frac);
     type1resHist->Fill(mean_type1_residual);
     type2resHist->Fill(mean_type2_residual);
 
     for (unsigned int ils = lumiInfoMapIterator->first.first; ils < lumiInfoMapIterator->first.second + 1; ils++) {
       if (ils > maxLS) {
         std::cout << "ils out of maxLS range!!" << std::endl;
         break;
       }
       Type1FracMon->setBinContent(ils, type1Frac);
       Type1FracMon->setBinError(ils, mean_type1_residual_unc);
       Type1ResMon->setBinContent(ils, mean_type1_residual);
       Type1ResMon->setBinError(ils, mean_type1_residual_unc);
       Type2ResMon->setBinContent(ils, mean_type2_residual);
       Type2ResMon->setBinError(ils, mean_type2_residual_unc);
     }
 
     type1FracGraph->SetPoint(iBlock, thisIOV + approxLumiBlockSize_ / 2.0, type1Frac);
     type1resGraph->SetPoint(iBlock, thisIOV + approxLumiBlockSize_ / 2.0, mean_type1_residual);
     type2resGraph->SetPoint(iBlock, thisIOV + approxLumiBlockSize_ / 2.0, mean_type2_residual);
     type1FracGraph->SetPointError(iBlock, approxLumiBlockSize_ / 2.0, mean_type1_residual_unc);
     type1resGraph->SetPointError(iBlock, approxLumiBlockSize_ / 2.0, mean_type1_residual_unc);
     type2resGraph->SetPointError(iBlock, approxLumiBlockSize_ / 2.0, mean_type2_residual_unc);
 
     edm::LogInfo("INFO")
         << "iBlock type1Frac mean_type1_residual mean_type2_residual mean_type1_residual_unc mean_type2_residual_unc "
         << iBlock << " " << type1Frac << " " << mean_type1_residual << " " << mean_type2_residual << " "
         << mean_type1_residual_unc << " " << mean_type2_residual_unc;
 
     type1Frac = 0.0;
     mean_type1_residual = 0.0;
     mean_type2_residual = 0.0;
     mean_type1_residual_unc = 0;
     mean_type2_residual_unc = 0;
 
     iBlock++;
 
     resetBlock();
   }
   histoFile->cd();
   type1FracHist->Write();
   type1resHist->Write();
   type2resHist->Write();
 
   delete type1FracHist;
   delete type1resHist;
   delete type2resHist;
 
   delete[] histname1;
   delete[] histname2;
   delete[] histname3;
 
   for (lumiInfoMapIterator = lumiInfoMap.begin(); (lumiInfoMapIterator != lumiInfoMap.end()); ++lumiInfoMapIterator) {
     delete lumiInfoMapIterator->second;
   }
   for (unsigned int lumiSection = 0; lumiSection < lumiSections.size(); lumiSection++) {
     thisLS = lumiSections[lumiSection];
     delete lumiInfoMapPerLS[thisLS];
   }
   lumiInfoMap.clear();
   lumiInfoMapPerLS.clear();
   lumiSections.clear();
   lumiInfoCounter.clear();
 }

void CorrPCCProducer::estimateType1Frac	(	std::vector< float >	uncorrPCCPerBX,
		float &	type1Frac
	)

private

Definition at line 206 of file CorrPCCProducer.cc.

References correctionTemplate_, evaluateCorrectionResiduals(), mps_fire::i, gen::k, mean_type1_residual, and LumiConstants::numBX.

Referenced by calculateCorrections().

                                                                                          {
   std::vector<float> corrected_tmp_;
   for (size_t i = 0; i < uncorrPCCPerBX.size(); i++) {
     corrected_tmp_.push_back(uncorrPCCPerBX.at(i));
   }
 
   //Apply initial type 1 correction
   for (size_t k = 0; k < LumiConstants::numBX - 1; k++) {
     float bin_k = corrected_tmp_.at(k);
     corrected_tmp_.at(k + 1) = corrected_tmp_.at(k + 1) - type1Frac * bin_k;
   }
 
   //Apply type 2 correction
   for (size_t i = 0; i < LumiConstants::numBX - 1; i++) {
     for (size_t j = i + 1; j < i + LumiConstants::numBX - 1; j++) {
       float bin_i = corrected_tmp_.at(i);
       if (j < LumiConstants::numBX) {
         corrected_tmp_.at(j) = corrected_tmp_.at(j) - bin_i * correctionTemplate_.at(j - i);
       } else {
         corrected_tmp_.at(j - LumiConstants::numBX) =
             corrected_tmp_.at(j - LumiConstants::numBX) - bin_i * correctionTemplate_.at(j - i);
       }
     }
   }
 
   //Apply additional iteration for type 1 correction
   evaluateCorrectionResiduals(corrected_tmp_);
   type1Frac += mean_type1_residual;
 }

void CorrPCCProducer::evaluateCorrectionResiduals ( std::vector< float > corrected_tmp_ )

private

Definition at line 237 of file CorrPCCProducer.cc.

References getMaximum(), histoFile, mean_type1_residual, mean_type1_residual_unc, mean_type2_residual, mean_type2_residual_unc, nTrain, LumiConstants::numBX, and electronIdCutBased_cfi::threshold.

Referenced by calculateCorrections(), and estimateType1Frac().

                                                                                  {
   float lumiMax = getMaximum(corrected_tmp_);
   float threshold = lumiMax * 0.2;
 
   mean_type1_residual = 0;
   mean_type2_residual = 0;
   nTrain = 0;
   float lumi = 0;
   std::vector<float> afterGlow;
   TH1F type1("type1", "", 1000, -0.5, 0.5);
   TH1F type2("type2", "", 1000, -0.5, 0.5);
   for (size_t ibx = 2; ibx < LumiConstants::numBX - 5; ibx++) {
     lumi = corrected_tmp_.at(ibx);
     afterGlow.clear();
     afterGlow.push_back(corrected_tmp_.at(ibx + 1));
     afterGlow.push_back(corrected_tmp_.at(ibx + 2));
 
     //Where type 1 and type 2 residuals are computed
     if (lumi > threshold && afterGlow[0] < threshold && afterGlow[1] < threshold) {
       for (int index = 3; index < 6; index++) {
         float thisAfterGlow = corrected_tmp_.at(ibx + index);
         if (thisAfterGlow < threshold) {
           afterGlow.push_back(thisAfterGlow);
         } else {
           break;
         }
       }
       float thisType1 = 0;
       float thisType2 = 0;
       if (afterGlow.size() > 1) {
         int nAfter = 0;
         for (unsigned int index = 1; index < afterGlow.size(); index++) {
           thisType2 += afterGlow[index];
           type2.Fill(afterGlow[index] / lumi);
           nAfter++;
         }
         thisType2 /= nAfter;
       }
       thisType1 = (afterGlow[0] - thisType2) / lumi;
 
       type1.Fill(thisType1);
       nTrain += 1;
     }
   }
 
   mean_type1_residual = type1.GetMean();  //Calculate the mean value of the type 1 residual
   mean_type2_residual = type2.GetMean();  //Calculate the mean value of the type 2 residual
   mean_type1_residual_unc = type1.GetMeanError();
   mean_type2_residual_unc = type2.GetMeanError();
 
   histoFile->cd();
   type1.Write();
   type2.Write();
 }

float CorrPCCProducer::getMaximum ( std::vector< float > lumi_vector )

private

Definition at line 194 of file CorrPCCProducer.cc.

References mps_fire::i.

Referenced by calculateCorrections(), and evaluateCorrectionResiduals().

                                                               {
   float max_lumi = 0;
   for (size_t i = 0; i < lumi_vector.size(); i++) {
     if (lumi_vector.at(i) > max_lumi)
       max_lumi = lumi_vector.at(i);
   }
   return max_lumi;
 }

void CorrPCCProducer::makeCorrectionTemplate ( )

private

Definition at line 186 of file CorrPCCProducer.cc.

References correctionTemplate_, JetChargeProducer_cfi::exp, LumiConstants::numBX, type2_a_, and type2_b_.

Referenced by CorrPCCProducer().

                                              {
   for (unsigned int bx = 1; bx < LumiConstants::numBX; bx++) {
     correctionTemplate_.at(bx) = type2_a_ * exp(-(float(bx) - 1) * type2_b_);
   }
 }

void CorrPCCProducer::resetBlock ( )

private

Definition at line 634 of file CorrPCCProducer.cc.

References correctionScaleFactors_, events_, LumiConstants::numBX, totalLumiByBX_, and totalLumiByBXAvg_.

Referenced by CorrPCCProducer(), and endRunProduce().

                                  {
   for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++) {
     totalLumiByBX_[bx] = 0;
     totalLumiByBXAvg_[bx] = 0;
     events_[bx] = 0;
     correctionScaleFactors_[bx] = 1.0;
   }
 }