Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

bsunanda:Run2-alca42 Update the trigger usage + provide tools for monitoring #14335

Merged
merged 2 commits into from Jun 3, 2016
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
26 changes: 17 additions & 9 deletions Calibration/HcalCalibAlgos/plugins/HcalIsoTrkAnalyzer.cc
Expand Up @@ -118,7 +118,7 @@ class HcalIsoTrkAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::
int t_Run, t_Event, t_ieta, t_goodPV, t_DataType;
double t_EventWeight, t_l1pt, t_l1eta, t_l1phi;
double t_l3pt, t_l3eta, t_l3phi, t_p, t_mindR1;
double t_mindR2, t_eMipDR, t_eHcal, t_hmaxNearP;
double t_mindR2, t_eMipDR, t_eHcal, t_eHcalDelta, t_hmaxNearP;
bool t_selectTk,t_qltyFlag,t_qltyMissFlag,t_qltyPVFlag;
std::vector<unsigned int> *t_DetIds;
std::vector<double> *t_HitEnergies, pbin;
Expand Down Expand Up @@ -483,6 +483,7 @@ void HcalIsoTrkAnalyzer::beginJob() {
tree->Branch("t_mindR2", &t_mindR2, "t_mindR2/D");
tree->Branch("t_eMipDR", &t_eMipDR, "t_eMipDR/D");
tree->Branch("t_eHcal", &t_eHcal, "t_eHcal/D");
tree->Branch("t_eHcalDelta", &t_eHcalDelta, "t_eHcalDelta/D");
tree->Branch("t_hmaxNearP", &t_hmaxNearP, "t_hmaxNearP/D");
tree->Branch("t_selectTk", &t_selectTk, "t_selectTk/O");
tree->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
Expand Down Expand Up @@ -631,13 +632,20 @@ int HcalIsoTrkAnalyzer::fillTree(std::vector< math::XYZTLorentzVector>& vecL1,
endcapRecHitsHandle, trkDetItr->pointHCAL,
trkDetItr->pointECAL, a_mipR_,
trkDetItr->directionECAL, nRH_eMipDR);
t_DetIds->clear(); t_HitEnergies->clear();
std::vector<DetId> ids;
t_eHcal = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
trkDetItr->pointECAL, a_coneR_,
trkDetItr->directionHCAL,nRecHits,
ids, *t_HitEnergies, useRaw_);
t_eHcal *= hcalScale_;
t_DetIds->clear(); t_HitEnergies->clear();
t_eHcalDelta = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
trkDetItr->pointECAL, a_charIsoR_,
trkDetItr->directionHCAL,nRecHits,
ids, *t_HitEnergies, useRaw_);
t_DetIds->clear(); t_HitEnergies->clear();
t_eHcal = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
trkDetItr->pointECAL, a_coneR_,
trkDetItr->directionHCAL,nRecHits,
ids, *t_HitEnergies, useRaw_);
t_eHcalDelta -= t_eHcal;
t_eHcalDelta *= hcalScale_;
t_eHcal *= hcalScale_;
for (unsigned int k=0; k<ids.size(); ++k) {
t_DetIds->push_back(ids[k].rawId());
}
Expand All @@ -650,8 +658,8 @@ int HcalIsoTrkAnalyzer::fillTree(std::vector< math::XYZTLorentzVector>& vecL1,
<< pTrack->phi() << "|" << t_p;
edm::LogInfo("HcalIsoTrack") << "e_MIP " << t_eMipDR
<< " Chg Isolation " << t_hmaxNearP
<< " eHcal" << t_eHcal << " ieta "
<< t_ieta << " Quality "
<< " eHcal" << t_eHcal << ":" << t_eHcalDelta
<< " ieta " << t_ieta << " Quality "
<< t_qltyMissFlag << ":"
<< t_qltyPVFlag << ":" << t_selectTk;
for (unsigned int lll=0;lll<t_DetIds->size();lll++) {
Expand Down
266 changes: 266 additions & 0 deletions Calibration/HcalCalibAlgos/test/CalibFitPlots.C
@@ -0,0 +1,266 @@
//////////////////////////////////////////////////////////////////////////////
// Usage:
// .L CalibFitPlots.C+g
// For standard set of histograms from CalibMonitor
// FitHistStandard(infile, outfile, prefix, mode, append, saveAll);
// For extended set of histograms from CalibMonitor
// FitHistExtended(infile, outfile, prefix, append);
//
// where:
// infile (std::string) = Name of the input ROOT file
// outfile (std::string) = Name of the output ROOT file
// prefix (std::string) = Prefix for the histogram names
// mode (int) = Flag to check which set of histograms to be
// done. It has the format lthdo where each of
// l, t,h,d,o can have a value 0 or 1 to select
// or deselect. l,t,h,d,o for momentum range
// 60-100, 30-40, all, 20-30, 40-60 Gev (11111)
// append (bool) = Open the output in Update/Recreate mode (True)
// saveAll (bool) = Flag to save intermediate plots (False)
//////////////////////////////////////////////////////////////////////////////

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TF1.h>
#include <TH1D.h>
#include <TProfile.h>
#include <TFitResult.h>
#include <TFitResultPtr.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TPaveStats.h>
#include <TPaveText.h>
#include <vector>
#include <string>
#include <iomanip>
#include <iostream>
#include <fstream>


std::pair<double,double> GetMean(TH1D* hist, double xmin, double xmax) {

double mean(0), rms(0), err(0), wt(0);
for (int i=1; i<=hist->GetNbinsX(); ++i) {
if (((hist->GetBinLowEdge(i)) >= xmin) ||
((hist->GetBinLowEdge(i)+hist->GetBinWidth(i)) <= xmax)) {
double cont = hist->GetBinContent(i);
double valu = hist->GetBinLowEdge(i)+0.5*+hist->GetBinWidth(i);
wt += cont;
mean += (valu*cont);
rms += (valu*valu*cont);
}
}
if (wt > 0) {
mean /= wt;
rms /= wt;
err = std::sqrt((rms-mean*mean)/wt);
}
return std::pair<double,double>(mean,err);
}

void FitHistStandard(std::string infile,std::string outfile,std::string prefix,
int mode=11111, bool append=true, bool saveAll=false) {

int iname[5] = {0, 1, 2, 3, 4};
int checkmode[5] = {10, 1000, 1, 10000, 100};
double xbins[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
double vbins[6] = {0.0, 7.0, 10.0, 13.0, 16.0, 50.0};
double dlbins[9] = {0.0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
std::string sname[4] = {"ratio","etaR", "dl1R","nvxR"};
std::string lname[4] = {"Z", "E", "L", "V"};
int numb[4] = {8, 8, 8, 5};
bool debug(true);

TFile *file = new TFile(infile.c_str());
std::vector<TH1D*> hists;
char name[100];
if (file != 0) {
for (int m1=0; m1<4; ++m1) {
for (int m2=0; m2<5; ++m2) {
sprintf (name, "%s%s%d0", prefix.c_str(), sname[m1].c_str(), iname[m2]);
TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
bool ok = ((hist0 != 0) && (hist0->GetEntries() > 25));
if ((mode/checkmode[m2])%10 > 0 && ok) {
TH1D* histo(0);
sprintf (name, "%s%s%d", prefix.c_str(), lname[m1].c_str(),iname[m2]);
if (m1 <= 1) histo = new TH1D(name, hist0->GetTitle(), numb[m1], xbins);
else if (m1 == 2) histo = new TH1D(name, hist0->GetTitle(), numb[m1], dlbins);
else histo = new TH1D(name, hist0->GetTitle(), numb[m1], vbins);
for (int j=0; j<=numb[m1]; ++j) {
sprintf (name, "%s%s%d%d", prefix.c_str(), sname[m1].c_str(), iname[m2], j);
TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
TH1D* hist = (TH1D*)hist1->Clone();
double value(0), error(0);
if (hist->GetEntries() > 0) {
value = hist->GetMean(); error = hist->GetRMS();
}
if (hist->GetEntries() > 4) {
double mean = hist->GetMean(), rms = hist->GetRMS();
double LowEdge = mean - 1.5*rms;
double HighEdge = mean + 2.0*rms;
if (LowEdge < 0.15) LowEdge = 0.15;
char option[20];
if (hist0->GetEntries() > 100) sprintf (option, "+QRS");
else sprintf (option, "+QRWLS");
double minvalue(0.30);
TFitResultPtr Fit = hist->Fit("gaus",option,"",LowEdge,HighEdge);
value = Fit->Value(1);
error = Fit->FitResult::Error(1);
std::pair<double,double> meaner = GetMean(hist,0.2,2.0);
if (debug) std::cout << "Fit " << value << ":" << error << ":"
<< hist->GetMeanError() << " Mean "
<< meaner.first << ":" << meaner.second;
if (value < minvalue || value > 2.0 || error > 0.5) {
value = meaner.first; error = meaner.second;
}
if (debug) std::cout << " Final " << value << ":" << error
<< std::endl;
}
if (j == 0) {
hists.push_back(hist);
} else {
if (saveAll) hists.push_back(hist);
histo->SetBinContent(j, value);
histo->SetBinError(j, error);
}
}
if (histo->GetEntries() > 2) {
int nbin = histo->GetNbinsX();
double LowEdge = histo->GetBinLowEdge(1);
double HighEdge= histo->GetBinLowEdge(nbin)+histo->GetBinWidth(nbin);
TFitResultPtr Fit = histo->Fit("pol0","+QRWLS","",LowEdge,HighEdge);
if (debug) std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- "
<< Fit->FitResult::Error(0) << std::endl;
histo->GetXaxis()->SetTitle("i#eta");
histo->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
histo->GetYaxis()->SetRangeUser(0.4,1.6);
}
hists.push_back(histo);
}
}
}
TFile* theFile(0);
if (append) {
theFile = new TFile(outfile.c_str(), "UPDATE");
} else {
theFile = new TFile(outfile.c_str(), "RECREATE");
}
theFile->cd();
for (unsigned int i=0; i<hists.size(); ++i) {
TH1D* hnew = (TH1D*)hists[i]->Clone();
hnew->Write();
}
theFile->Close();
file->Close();
}
}

void FitHistExtended(std::string infile,std::string outfile,std::string prefix,
bool append=true) {
double xbins[43]= {-21.5,-20.5,-19.5,-18.5,-17.5,-16.5,-15.5,-14.5,-13.5,
-12.5,-11.5,-10.5,-9.5,-8.5,-7.5,-6.5,-5.5,-4.5,-3.5,
-2.5,-1.5,0.0,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,
11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5};
std::string sname("ratio"), lname("Z");
int iname(2), numb(42);
bool debug(true);

TFile *file = new TFile(infile.c_str());
std::vector<TH1D*> hists;
char name[100];
if (file != 0) {
sprintf (name, "%s%s%d0", prefix.c_str(), sname.c_str(), iname);
TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
bool ok = (hist0 != 0);
if (ok) {
sprintf (name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
TH1D* histo = new TH1D(name, hist0->GetTitle(), numb, xbins);

int nbin = hist0->GetNbinsX();
if (hist0->GetEntries() > 10) {
double mean = hist0->GetMean(), rms = hist0->GetRMS();
double LowEdge = mean - 1.5*rms;
double HighEdge = mean + 2.0*rms;
if (LowEdge < 0.15) LowEdge = 0.15;
char option[20];
if (hist0->GetEntries() > 100) {
sprintf (option, "+QRS");
} else {
sprintf (option, "+QRWLS");
HighEdge= mean+1.5*rms;
}
TFitResultPtr Fit = hist0->Fit("gaus",option,"",LowEdge,HighEdge);
std::pair<double,double> meaner = GetMean(hist0,0.2,2.0);
if (debug) std::cout << "Fit " << Fit->Value(1) << ":"
<< Fit->FitResult::Error(1) << ":"
<< hist0->GetMeanError() << " Mean "
<< meaner.first << ":" << meaner.second;
}
for (int j=0; j<=numb; ++j) {
sprintf (name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j);
TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
TH1D* hist = (TH1D*)hist1->Clone();
double value(0), error(0), total(0);
if (hist->GetEntries() > 0) {
value = hist->GetMean(); error = hist->GetRMS();
for (int i=1; i<=nbin; ++i) total += hist->GetBinContent(i);
}
if (total > 4) {
double mean = hist->GetMean(), rms = hist->GetRMS();
double LowEdge = mean - 1.5*rms;
double HighEdge = mean + 2.0*rms;
if (LowEdge < 0.15) LowEdge = 0.15;
char option[20];
if (total > 100) {
sprintf (option, "+QRS");
} else {
sprintf (option, "+QRWLS");
HighEdge= mean+1.5*rms;
}
double minvalue(0.30);
TFitResultPtr Fit = hist->Fit("gaus",option,"",LowEdge,HighEdge);
value = Fit->Value(1);
error = Fit->FitResult::Error(1);
std::pair<double,double> meaner = GetMean(hist,0.2,2.0);
if (debug) std::cout << "Fit " << value << ":" << error << ":"
<< hist->GetMeanError() << " Mean "
<< meaner.first << ":" << meaner.second;
if (value < minvalue || value > 2.0 || error > 0.5) {
value = meaner.first; error = meaner.second;
}
if (debug) std::cout << " Final " << value << ":" << error<<std::endl;
}
hists.push_back(hist);
if (j != 0) {
histo->SetBinContent(j, value);
histo->SetBinError(j, error);
}
}
if (histo->GetEntries() > 2) {
TFitResultPtr Fit = histo->Fit("pol0","+QRWLS","",xbins[0],xbins[numb]);
if (debug) std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- "
<< Fit->FitResult::Error(0) << std::endl;
histo->GetXaxis()->SetTitle("i#eta");
histo->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
histo->GetYaxis()->SetRangeUser(0.4,1.6);
}
hists.push_back(histo);
}
TFile* theFile(0);
if (append) {
theFile = new TFile(outfile.c_str(), "UPDATE");
} else {
theFile = new TFile(outfile.c_str(), "RECREATE");
}

theFile->cd();
for (unsigned int i=0; i<hists.size(); ++i) {
TH1D* hnew = (TH1D*)hists[i]->Clone();
hnew->Write();
}
theFile->Close();
file->Close();
}
}