Skip to content

Commit

Permalink
Merge pull request #41781 from mmusich/migrate_DiMuonMassBiasClient_t…
Browse files Browse the repository at this point in the history 
…o_Profiles_13_1_X

[13.1.X] Migrate `DiMuonMassBiasClient` mass bias histograms to profiles
  • Loading branch information
@cmsbuild
cmsbuild authored May 28, 2023
2 parents 92b074f + 37ecc1b commit 2fcc6f8
Show file tree
Hide file tree
Showing 3 changed files with 167 additions and 35 deletions.
7 changes: 7 additions & 0 deletions DQMOffline/Alignment/interface/DiMuonMassBiasClient.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,9 +81,12 @@ class DiMuonMassBiasClient : public DQMEDHarvester {
void bookMEs(DQMStore::IBooker& ibooker);
void getMEsToHarvest(DQMStore::IGetter& igetter);
diMuonMassBias::fitOutputs fitLineShape(TH1* hist, const bool& fitBackground = false) const;
void fitAndFillProfile(std::pair<std::string, MonitorElement*> toHarvest, DQMStore::IBooker& iBooker);
void fitAndFillHisto(std::pair<std::string, MonitorElement*> toHarvest, DQMStore::IBooker& iBooker);

// data members
const std::string TopFolder_;
const bool useTH1s_;
const bool fitBackground_;
const bool useRooCBShape_;
const bool useRooCMSShape_;
Expand All @@ -97,6 +100,10 @@ class DiMuonMassBiasClient : public DQMEDHarvester {
std::vector<std::string> MEtoHarvest_;

// the histograms to be filled
std::map<std::string, MonitorElement*> meanHistos_;
std::map<std::string, MonitorElement*> widthHistos_;

// the profiles to be filled
std::map<std::string, MonitorElement*> meanProfiles_;
std::map<std::string, MonitorElement*> widthProfiles_;

Expand Down
3 changes: 2 additions & 1 deletion DQMOffline/Alignment/python/ALCARECOTkAlDQM_cff.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,8 @@

ALCARECOTkAlDiMuonMassBiasDQM = DQMOffline.Alignment.DiMuonMassBiasMonitor_cfi.DiMuonMassBiasMonitor.clone(
muonTracks = 'ALCARECO'+__trackCollName,
FolderName = "AlCaReco/"+__selectionName
FolderName = "AlCaReco/"+__selectionName,
DiMuMassConfig = dict(maxDeltaEta = 3.5)
)

ALCARECOTkAlDiMuonAndVertexDQM = cms.Sequence(ALCARECOTkAlDiMuonAndVertexTkAlDQM + ALCARECOTkAlDiMuonAndVertexVtxDQM + ALCARECOTkAlDiMuonMassBiasDQM)
Expand Down
192 changes: 158 additions & 34 deletions DQMOffline/Alignment/src/DiMuonMassBiasClient.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
//-----------------------------------------------------------------------------------
DiMuonMassBiasClient::DiMuonMassBiasClient(edm::ParameterSet const& iConfig)
: TopFolder_(iConfig.getParameter<std::string>("FolderName")),
useTH1s_(iConfig.getParameter<bool>("useTH1s")),
fitBackground_(iConfig.getParameter<bool>("fitBackground")),
useRooCBShape_(iConfig.getParameter<bool>("useRooCBShape")),
useRooCMSShape_(iConfig.getParameter<bool>("useRooCMSShape")),
Expand Down Expand Up @@ -85,12 +86,12 @@ void DiMuonMassBiasClient::bookMEs(DQMStore::IBooker& iBooker)
const auto& xmax = ME->getAxisMax(1);

MonitorElement* meanToBook =
iBooker.book1D(("Mean" + key), (title + ";" + xtitle + ";" + ytitle), nxbins, xmin, xmax);
meanProfiles_.insert({key, meanToBook});
iBooker.book1D(("Mean" + key), (title + ";#LT M_{#mu^{-}#mu^{+}} #GT [GeV];" + ytitle), nxbins, xmin, xmax);
meanHistos_.insert({key, meanToBook});

MonitorElement* sigmaToBook =
iBooker.book1D(("Sigma" + key), (title + ";" + xtitle + ";" + "#sigma of " + ytitle), nxbins, xmin, xmax);
widthProfiles_.insert({key, sigmaToBook});
widthHistos_.insert({key, sigmaToBook});
}
}

Expand All @@ -113,46 +114,168 @@ void DiMuonMassBiasClient::getMEsToHarvest(DQMStore::IGetter& iGetter)
}
}

//-----------------------------------------------------------------------------------
void DiMuonMassBiasClient::fitAndFillProfile(std::pair<std::string, MonitorElement*> toHarvest,
DQMStore::IBooker& iBooker)
//-----------------------------------------------------------------------------------
{
const auto& key = toHarvest.first;
const auto& ME = toHarvest.second;

if (debugMode_)
edm::LogPrint("DiMuonMassBiasClient") << "dealing with key: " << key << std::endl;

if (ME == nullptr) {
edm::LogError("DiMuonMassBiasClient") << "could not find MonitorElement for key: " << key << std::endl;
return;
}

const auto& title = ME->getTitle();
const auto& xtitle = ME->getAxisTitle(1);
const auto& ytitle = ME->getAxisTitle(2);

const auto& nxbins = ME->getNbinsX();
const auto& xmin = ME->getAxisMin(1);
const auto& xmax = ME->getAxisMax(1);

TProfile* p_mean = new TProfile(("Mean" + key).c_str(),
(title + ";" + xtitle + ";#LT M_{#mu^{-}#mu^{+}} #GT [GeV]").c_str(),
nxbins,
xmin,
xmax,
"g");

TProfile* p_width = new TProfile(
("Sigma" + key).c_str(), (title + ";" + xtitle + ";#sigma of " + ytitle).c_str(), nxbins, xmin, xmax, "g");

p_mean->Sumw2();
p_width->Sumw2();

TH2F* bareHisto = ME->getTH2F();
for (int bin = 1; bin <= nxbins; bin++) {
const auto& xaxis = bareHisto->GetXaxis();
const auto& low_edge = xaxis->GetBinLowEdge(bin);
const auto& high_edge = xaxis->GetBinUpEdge(bin);

if (debugMode_)
edm::LogPrint("DiMuonMassBiasClient") << "dealing with bin: " << bin << " range: (" << std::setprecision(2)
<< low_edge << "," << std::setprecision(2) << high_edge << ")";

TH1D* Proj = bareHisto->ProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin);
Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin));

diMuonMassBias::fitOutputs results = fitLineShape(Proj);

if (results.isInvalid()) {
edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl;
continue;
}

// fill the mean profiles
const Measurement1D& bias = results.getBias();

// ============================================= DISCLAIMER ================================================
// N.B. this is sort of a hack in order to fill arbitrarily both central values and error bars of a TProfile.
// Choosing the option "g" in the constructor the bin error will be 1/sqrt(W(j)), where W(j) is the sum of weights.
// Filling the sum of weights with the 1 / err^2, the bin error automatically becomes "err".
// In order to avoid the central value to be shifted, that's divided by 1 / err^2 as well.
// For more information, please consult the https://root.cern.ch/doc/master/classTProfile.html

p_mean->SetBinContent(bin, bias.value() / (bias.error() * bias.error()));
p_mean->SetBinEntries(bin, 1. / (bias.error() * bias.error()));

if (debugMode_)
LogDebug("DiMuonBassBiasClient") << " Bin: " << bin << " value: " << bias.value() << " from profile ( "
<< p_mean->GetBinContent(bin) << ") - error: " << bias.error()
<< " from profile ( " << p_mean->GetBinError(bin) << " )";

// fill the width profiles
const Measurement1D& width = results.getWidth();

// see discussion above
p_width->SetBinContent(bin, width.value() / (width.error() * width.error()));
p_width->SetBinEntries(bin, 1. / (width.error() * width.error()));
}

// now book the profiles
iBooker.setCurrentFolder(TopFolder_ + "/DiMuonMassBiasMonitor/MassBias/Profiles");
MonitorElement* meanToBook = iBooker.bookProfile(p_mean->GetName(), p_mean);
meanProfiles_.insert({key, meanToBook});

MonitorElement* sigmaToBook = iBooker.bookProfile(p_width->GetName(), p_width);
widthProfiles_.insert({key, sigmaToBook});

delete p_mean;
delete p_width;
}

//-----------------------------------------------------------------------------------
void DiMuonMassBiasClient::fitAndFillHisto(std::pair<std::string, MonitorElement*> toHarvest,
DQMStore::IBooker& iBooker)
//-----------------------------------------------------------------------------------
{
const auto& key = toHarvest.first;
const auto& ME = toHarvest.second;

if (debugMode_)
edm::LogPrint("DiMuonMassBiasClient") << "dealing with key: " << key << std::endl;

if (ME == nullptr) {
edm::LogError("DiMuonMassBiasClient") << "could not find MonitorElement for key: " << key << std::endl;
return;
}

TH2F* bareHisto = ME->getTH2F();
for (int bin = 1; bin <= ME->getNbinsX(); bin++) {
const auto& xaxis = bareHisto->GetXaxis();
const auto& low_edge = xaxis->GetBinLowEdge(bin);
const auto& high_edge = xaxis->GetBinUpEdge(bin);

if (debugMode_)
edm::LogPrint("DiMuonMassBiasClient") << "dealing with bin: " << bin << " range: (" << std::setprecision(2)
<< low_edge << "," << std::setprecision(2) << high_edge << ")";
TH1D* Proj = bareHisto->ProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin);
Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin));

diMuonMassBias::fitOutputs results = fitLineShape(Proj);

if (results.isInvalid()) {
edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl;
continue;
}

// fill the mean profiles
const Measurement1D& bias = results.getBias();
meanHistos_[key]->setBinContent(bin, bias.value());
meanHistos_[key]->setBinError(bin, bias.error());

// fill the width profiles
const Measurement1D& width = results.getWidth();
widthHistos_[key]->setBinContent(bin, width.value());
widthHistos_[key]->setBinError(bin, width.error());
}
}

//-----------------------------------------------------------------------------------
void DiMuonMassBiasClient::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter)
//-----------------------------------------------------------------------------------
{
edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::endLuminosityBlock";

getMEsToHarvest(igetter);
bookMEs(ibooker);

for (const auto& [key, ME] : harvestTargets_) {
if (debugMode_)
edm::LogPrint("DiMuonMassBiasClient") << "dealing with key: " << key << std::endl;
TH2F* bareHisto = ME->getTH2F();
for (int bin = 1; bin <= ME->getNbinsX(); bin++) {
const auto& xaxis = bareHisto->GetXaxis();
const auto& low_edge = xaxis->GetBinLowEdge(bin);
const auto& high_edge = xaxis->GetBinUpEdge(bin);

if (debugMode_)
edm::LogPrint("DiMuonMassBiasClient") << "dealing with bin: " << bin << " range: (" << std::setprecision(2)
<< low_edge << "," << std::setprecision(2) << high_edge << ")";
TH1D* Proj = bareHisto->ProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin);
Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin));

diMuonMassBias::fitOutputs results = fitLineShape(Proj);

if (results.isInvalid()) {
edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl;
continue;
}

// fill the mean profiles
const Measurement1D& bias = results.getBias();
meanProfiles_[key]->setBinContent(bin, bias.value());
meanProfiles_[key]->setBinError(bin, bias.error());
// book the histograms upfront
if (useTH1s_) {
bookMEs(ibooker);
}

// fill the width profiles
const Measurement1D& width = results.getWidth();
widthProfiles_[key]->setBinContent(bin, width.value());
widthProfiles_[key]->setBinError(bin, width.error());
for (const auto& element : harvestTargets_) {
if (!useTH1s_) {
// if using profiles
this->fitAndFillProfile(element, ibooker);
} else {
// if using histograms
this->fitAndFillHisto(element, ibooker);
}
}
}
Expand Down Expand Up @@ -298,6 +421,7 @@ void DiMuonMassBiasClient::fillDescriptions(edm::ConfigurationDescriptions& desc
{
edm::ParameterSetDescription desc;
desc.add<std::string>("FolderName", "DiMuonMassBiasMonitor");
desc.add<bool>("useTH1s", false);
desc.add<bool>("fitBackground", false);
desc.add<bool>("useRooCMSShape", false);
desc.add<bool>("useRooCBShape", false);
Expand Down

0 comments on commit 2fcc6f8

Please sign in to comment.