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

new function of clipping when generating L2 corrections, better fit used plotting L1 offsets. #14

Open
wants to merge 12 commits into
base: master
Choose a base branch
from
Open
199 changes: 101 additions & 98 deletions JetAnalyzers/bin/compareJEC.C

Large diffs are not rendered by default.

208 changes: 104 additions & 104 deletions JetAnalyzers/bin/jet_synchplot_x.cc

Large diffs are not rendered by default.

22 changes: 12 additions & 10 deletions JetAnalyzers/bin/jet_synchtest_x.cc
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ class MatchEventsAndJets {
int NBinsNpvRhoNpu;
int npvRhoNpuBinWidth;
int iIT;
int inpv;
int inpv;
int inpv_low;
int inpv_high;
int irho;
Expand Down Expand Up @@ -339,9 +339,10 @@ void MatchEventsAndJets::WriteMatchedEventsMaps(map<evtid, pair<ull,ull>, evtid>
outputFilename = "matchedEventsMaps_"+algo1+".root";
outputFilename = outputPath+outputFilename;
string name = "mapTree";
cout << "\tWriting " << name+(noPU ? "NoPU" : "PU") << " to " << outputFilename << " ... " << flush;
cout << "\tWriting " << name+(noPU ? "NoPU" : "PU") << " to " << outputFilename << " ... " << flush;
string option = (noPU ? "UPDATE" : "RECREATE");
TFile* mapFile = TFile::Open(outputFilename.c_str(),option.c_str());
mapFile->cd();
mapFile->WriteObject(&mapTree,(name+(noPU ? "NoPU" : "PU")).c_str());
mapFile->Write();
mapFile->Close();
Expand Down Expand Up @@ -440,7 +441,7 @@ void MatchEventsAndJets::DeclareHistograms(bool reduceHistograms) {
histograms["g_GenWeight"] = new TH1D("g_GenWeight", "g_GenWeight;log_{10}(GenWeight);Events", 1000,-48,2);
histograms["g_pThatWeight"] = new TH1D("g_pThatWeight;log_{10}(pThatWeight);Events","g_pThatWeight", 1000,-48,2);
histograms["g_weight"] = new TH1D("g_weight","g_weight;log_{10}(EvtWeight);Events", 1000,-48,2);
histograms["g_pthat"] = new TH1D("g_pthat","g_pthat;#hat{p}_{T}^{PU};Events",(int)vpt[NPtBins]/10.0,vpt[0],vpt[NPtBins]);
histograms["g_pthat"] = new TH1D("g_pthat","g_pthat;#hat{p}_{T}^{PU};Events",(int)vpt[NPtBins]/10.0,vpt[0],vpt[NPtBins]);
if(!reduceHistograms) {
histograms["g_nj"] = new TH2D("g_nj","g_nj",30,0,30,30,0,30);
histograms["g_npv"] = new TH2D("g_npv","g_npv",50,0,50,50,0,50);
Expand Down Expand Up @@ -795,6 +796,7 @@ void MatchEventsAndJets::LoopOverEvents(bool verbose, bool reduceHistograms, str
for (IT::const_iterator it = mapTreePU.begin(); it != mapTreePU.end(); ++it) {

if (iftest && nevs >= maxEvts) return;
// if (nevs >= 10000000) return;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove this line

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done


//if (nevs%10000==0) cout << "\t"<<nevs << endl;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove this line too.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

loadbar2(nevs+1,nentries,50,"\t\t");
Expand Down Expand Up @@ -830,7 +832,7 @@ void MatchEventsAndJets::LoopOverEvents(bool verbose, bool reduceHistograms, str
else {
jetMapIndex++;
ReadJetMap(jetMapIndex,readJetMap);
}
}

if(FillHistograms(reduceHistograms)) nevs++;

Expand Down Expand Up @@ -908,7 +910,7 @@ void MatchEventsAndJets::FillRecToRecThroughGenMap() {
}
}
if(j1 >= 0 && j2 >= 0 && j1 < tpu->nref && j2 < tnopu->nref &&
tpu->refdrjt->at(j1) < maxDeltaR && tnopu->refdrjt->at(j2) < maxDeltaR &&
tpu->refdrjt->at(j1) < maxDeltaR && tnopu->refdrjt->at(j2) < maxDeltaR &&
fabs(tpu->refpt->at(j1) - tnopu->refpt->at(j2))<0.0001) {
jetMap[j1] = j2;
}
Expand Down Expand Up @@ -1186,8 +1188,8 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) {
double avg_offset = 0;
double avg_offset_det[NDetectorNames] = {0,0,0,0};
double njet_det[NDetectorNames] = {0,0,0,0};
// MATCHING HISTOS.

// MATCHING HISTOS.
// Loop over matched jets
int jpu = -1;
int jnopu = -1;
Expand Down Expand Up @@ -1356,7 +1358,7 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) {
dynamic_cast<TH2D*>(histograms[hname])->Fill(tpu->refpt->at(jpu),respNopu,weight);
hname = Form("p_offAfterOoffBeforeVsrefpt_%s_npv%i_%i",detectorAbbreviation.Data(),inpv_low,inpv_high);
dynamic_cast<TH2D*>(histograms[hname])->Fill(tpu->refpt->at(jpu),offset/offset_raw,weight);

//RHO
hname = Form("p_resVsrefpt_%s_rho%i_%i",detectorAbbreviation.Data(),irho_low,irho_high);
dynamic_cast<TH2D*>(histograms[hname])->Fill(tpu->refpt->at(jpu),resp,weight);
Expand All @@ -1372,7 +1374,7 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) {
dynamic_cast<TH2D*>(histograms[hname])->Fill(tpu->refpt->at(jpu),respNopu,weight);
hname = Form("p_offAfterOoffBeforeVsrefpt_%s_rho%i_%i",detectorAbbreviation.Data(),irho_low,irho_high);
dynamic_cast<TH2D*>(histograms[hname])->Fill(tpu->refpt->at(jpu),offset/offset_raw,weight);

//OTHER
hname = Form("p_resVsnpu_%s_pt%.1f_%.1f",detectorAbbreviation.Data(),
vpt[JetInfo::getBinIndex(tpu->refpt->at(jpu),vpt,NPtBins)],vpt[JetInfo::getBinIndex(tpu->refpt->at(jpu),vpt,NPtBins)+1]);
Expand All @@ -1389,7 +1391,7 @@ bool MatchEventsAndJets::FillHistograms(bool reduceHistograms) {
avg_offset_det[idet]+=offset;
njet_det[idet]+=1.;

} // for matched jets
} // for matched jets

avg_offset /= jetMap.size();
for (int det=0;det<NDetectorNames;det++) {
Expand Down
6 changes: 3 additions & 3 deletions JetAnalyzers/test/run_JRA_cfg.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
#!
# Conditions source options: GT, SQLite, DB
conditionsSource = "GT"
era = "Spring16_25nsV1_MC"
era = "Fall17_25nsV1_MC"
doProducer = False
process = cms.Process("JRA")
multithread = False
Expand Down Expand Up @@ -48,13 +48,13 @@
#! CONDITIONS (DELIVERING JEC BY DEFAULT!)
#!
process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff")
process.GlobalTag.globaltag = cms.string('80X_mcRun2_asymptotic_v5_2016PixDynIneff')
process.GlobalTag.globaltag = cms.string('94X_mc2017_realistic_v10')

if conditionsSource != "GT":
if conditionsSource == "DB":
conditionsConnect = cms.string("frontier://FrontierPrep/CMS_COND_PHYSICSTOOLS")
elif conditionsSource == "SQLite":
conditionsConnect = cms.string('sqlite_file:'+era+'.db')
conditionsConnect = cms.string('sqlite_file:'+era+'.db')

from CondCore.DBCommon.CondDBSetup_cfi import *
process.jec = cms.ESSource("PoolDBESSource",CondDBSetup,
Expand Down
2 changes: 0 additions & 2 deletions JetUtilities/bin/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -12,5 +12,3 @@
</bin>
<bin name="jet_add_histos_x" file="jet_add_histos_x.cc">
</bin>
<bin name="test_jet_corrector_x" file="test_jet_corrector_x.cc">
</bin>
9 changes: 5 additions & 4 deletions JetUtilities/interface/ClosureMaker.hh
Original file line number Diff line number Diff line change
Expand Up @@ -83,27 +83,28 @@ public:
void checkResponse();
pair<double,double> determineCanvasRange(double xmin, double xmax);
void makeCanvases();
void makeMergedCanvas();
void makeMergedCanvas(bool finemerge);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any way to fix the spacing here? I know it's not your fault, but this would be a good opportunity to standardize the use of spaces and tabs.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am converting all tabs to spaces when I modify a file now.

void writeToFile();
void makeClosure(const VARIABLES::Variable ivar = VARIABLES::refpt);

private:
bool objects_loaded, draw_guidelines;
double CMEnergy, nsigma;
double CMEnergy, nsigma;
TString path, filename, outputDir, outputFilename, flavor, alg, histMet;
vector<TString> algs, outputFormat;
JetInfo *ji;
TFile *ifile, *ofile;
ObjectLoader<TH2F> hl;
vector<TH1D*> h;
vector<TF1*> func;
vector<TH1D*> h;
vector<TF1*> func;
vector<TH1F*> hClosure;
TF1 *line, *linePlus, *lineMinus;
vector<pair<TCanvas*,TLegend*> > canvases_legends;
vector<TPaveText*> pave;
VARIABLES::Variable var;
TDirectoryFile *odir;
HistUtil::HistogramMetric histogramMetric;
int statTh;
};

#endif
11 changes: 8 additions & 3 deletions JetUtilities/interface/HistogramUtilities.hh
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include "TMath.h"
#include "TH2.h"
#include "TF1.h"
#include "TSpectrum.h"

//C++ libraries
#include <iostream>
Expand Down Expand Up @@ -48,12 +49,12 @@ namespace HistUtil{
// (2) The mean (mu_f or mpv) and RMS (RMS_f or sigma_f) of a fit to the histogram (i.e. a Gaussian fit)
// (3) The median of the histogram
enum HistogramMetric {none, mu_h, RMS_h, mu_f, mpv, RMS_f, sigma_f, median};
static const unsigned int nHistogramMetric = 8;
static const unsigned int nHistogramMetric = 8;

enum AxisDirection {X, Y, Z};
static const unsigned int nAxisDirections = 3;

// A routine that returns the string given the HistogramMetric
// A routine that returns the string given the HistogramMetric
string getHistogramMetricString(HistogramMetric);

// A routine that returns the HistogramMetric type given the string
Expand All @@ -64,7 +65,7 @@ namespace HistUtil{

// A routine that returns a given HistogramMetric and its error (if any)
pair<double,double> getHistogramMetric1D(HistogramMetric, TH1*, double fallback_threshold = 0.05, bool verbose = false);
pair<double,double> getHistogramMetric1D(string, TH1*, double fallback_threshold = 0.05, bool verbose = false);
pair<double,double> getHistogramMetric1D(string, TH1*, double fallback_threshold = 0.05, bool verbose = false);

// A routine that returns the median of a given histogram and an error equal to the width of the bin that contains the histogram
pair<double,double> getHistogramMedian1D(TH1*, bool debug = false);
Expand All @@ -91,6 +92,10 @@ namespace HistUtil{
// (i.e. if the AxisDirection=Y, then the metric will be plotted along X)
vector<pair<double,double> > getHistogramMetric2D(HistogramMetric, TH2*, AxisDirection, const vector<string>&);
vector<pair<double,double> > getHistogramMetric2D(string, TH2*, AxisDirection, const vector<string>&);

void adjust_fitrange(TH1* h,double& min,double& max);
int number_filled_bins(TH1* h, double min, double max);
TF1* fit_gaussian(TH1*& hrsp, const double nsigma, const double jtptmin, const int niter, const int verbose);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you're going to add fit_gaussian here why not also add fit_dscb?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because I am not sure what that is. And I have never got a chance to use this fit_dscb() yet.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a fit for a double sided crystal ball function. It performs the same type of operation as the fit_gaussian function, but for a different formula.

}

#endif
3 changes: 2 additions & 1 deletion JetUtilities/interface/L2Creator.hh
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ private:
TString output, outputDir, l2calofit, l2pffit;
vector<string> formats, algs;
bool l2l3, delphes;
int maxFitIter;
int maxFitIter, statTh;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It wasn't immediately clear to me from the name what this variable was. Can you call it "statThreshold"?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

HistUtil::HistogramMetric histogramMetric;
TFile* ofile;
TFile* ifile;
Expand All @@ -116,6 +116,7 @@ private:
vector<TGraphErrors*> vabscor_eta;
vector<TGraph*> vrelcor_eta;
vector<PiecewiseSpline*> vabscor_eta_spline;
float ptclip;
};

#endif
38 changes: 19 additions & 19 deletions JetUtilities/src/ClosureMaker.cc
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@ ClosureMaker::ClosureMaker(CommandLine& cl) {
histMet = cl.getValue<TString> ("histMet", "mu_h");
histogramMetric = HistUtil::getHistogramMetricType(string(histMet));
bool help = cl.getValue<bool> ("help", false);
statTh = cl.getValue<int> ("statTh", 4);

if (help) {cl.print(); return;}
if (!cl.partialCheck()) return;
Expand Down Expand Up @@ -88,7 +89,7 @@ void ClosureMaker::openInputFile() {
ifile = TFile::Open(path+filename+"_"+flavor+".root","READ");
else
ifile = TFile::Open(path+filename+".root","READ");

if(ifile == 0) {
cout << "ERROR::ClosureMaker::openInputFiles() Could not open the file " << path << endl;
std::terminate();
Expand All @@ -102,7 +103,7 @@ void ClosureMaker::getHistograms(TDirectoryFile* idir) {
}

if(var == VARIABLES::refpt || var == VARIABLES::ptclpt) {
hl.load_objects(idir,"RelRspVsRefPt:JetEta");
hl.load_objects(idir,"RelRspVsRefPt:JetEta");
objects_loaded = true;
if(hl.nobjects()!=NDetectorNames) {
cout << "One or more of the histogram pointers from file " << path << " is NULL." << endl;
Expand All @@ -125,7 +126,7 @@ void ClosureMaker::getHistograms(TDirectoryFile* idir) {
void ClosureMaker::openOutputFile() {
//
// Open/create the output directory and file
//
//
TString ofname = Form("%s/ClosureVs%s.root",outputDir.Data(),getVariableTitleString(var).c_str());
if(!flavor.IsNull()) ofname = Form("%s/ClosureVs%s_%s.root",outputDir.Data(),
getVariableTitleString(var).c_str(),flavor.Data());
Expand Down Expand Up @@ -209,7 +210,7 @@ void ClosureMaker::loopOverAlgorithms() {
if (strcmp(dirKey->GetClassName(),"TDirectoryFile")!=0) continue;
TDirectoryFile* idir = (TDirectoryFile*)dirKey->ReadObj();
alg = idir->GetName(); if (!JetInfo::contains(algs,alg)) continue;

algIndex++;
cout << alg << " ... " << endl;

Expand Down Expand Up @@ -258,7 +259,7 @@ void ClosureMaker::loopOverAlgorithms() {
//
makeCanvases();
if(var == VARIABLES::refpt || var == VARIABLES::jtpt || var == VARIABLES::ptclpt) {
makeMergedCanvas();
makeMergedCanvas(false);
}

writeToFile();
Expand Down Expand Up @@ -343,7 +344,7 @@ void ClosureMaker::loopOverBins(TH2F* hvar, unsigned int iVarBin) {
varBins[ibin].c_str(),varBins[ibin+1].c_str());
h.push_back(hvar->ProjectionY(name,ibin+1,ibin+1,"e"));

if (h.back()->GetEntries()>4) {
if (h.back()->GetEntries()> statTh) {
if(histogramMetric==HistUtil::mu_f || histogramMetric==HistUtil::mpv) {
int nbins = 50;//100;
TSpectrum *spec = new TSpectrum(10);
Expand Down Expand Up @@ -385,10 +386,10 @@ void ClosureMaker::loopOverBins(TH2F* hvar, unsigned int iVarBin) {
continue;
}
}
else if(h.back()->GetEntries()<=4 && h.back()->GetEntries()>1) {
hClosure.back()->SetBinContent(ibin+1,h.back()->GetMean());
hClosure.back()->SetBinError(ibin+1,h.back()->GetMeanError());
}
// else if(h.back()->GetEntries()<=4 && h.back()->GetEntries()>1) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why remove this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There was a point in the closure whose response plot had exactly 4 entries, and it was not removed by statThreshold option. So I removed this line

Copy link
Member

@aperloff aperloff Apr 9, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well yeah, that's the logic set up there. What do you want to have happen?

  1. If there are more than statThreshold entries and things proceed as normal
  2. The part you commented out. If there is less than statThreshold entries, but it's not empty, do you want some fallback behavior? That was what these lines were for. If you don't want that behavior don't just comment these lines out; remove them entirely.
  3. else skip this point

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If it does not exceed the statThreshold, I think it'd be good to skip the point to show that the correction is performing well. So I'll skip this point

// hClosure.back()->SetBinContent(ibin+1,h.back()->GetMean());
// hClosure.back()->SetBinError(ibin+1,h.back()->GetMeanError());
// }
else {
continue;
}
Expand Down Expand Up @@ -421,7 +422,7 @@ void ClosureMaker::checkResponse() {
cout << "\tWARNING Closure for " << hClosure[ih]->GetName() << " at "
<< hClosure[ih]->GetBinLowEdge(ibin) << " < "
<< getVariableTitleString(var) << " < " << hClosure[ih]->GetBinLowEdge(ibin+1)
<<" has relresp as low as " << binCont << " +/- " << binErr << endl;
<<" has relresp as low as " << binCont << " +/- " << binErr << endl;
}
}
}
Expand Down Expand Up @@ -491,7 +492,7 @@ void ClosureMaker::makeCanvases() {
tdrLeg(0.58,0.16,0.9,0.4)));
if((var == VARIABLES::refpt || var == VARIABLES::jtpt || var == VARIABLES::ptclpt) && ih<3)
canvases_legends.back().first->GetPad(0)->SetLogx();

//
// Format and draw the pave
//
Expand Down Expand Up @@ -533,10 +534,10 @@ void ClosureMaker::makeCanvases() {
}

//______________________________________________________________________________
void ClosureMaker::makeMergedCanvas() {
void ClosureMaker::makeMergedCanvas(bool finemerge = false) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't love this finemerge variable. To me it just seems like a way to force the canvas range to be narrower. Can't you just change the settings in the function "determineCanvasRange"? Make that function smarter rather than coming up with a way to circumvent it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, it is forcing the canvas to be narrowed to a range that high pT range should locate in. And also this is the plot people always want to see to make sure high pT region is fine. Same canvas range also makes it easier to compare between iterations.
There is such a determineCanvasRange function if I did not remember wrong. But when the closure it not good, it will try to incorporate all points in the range. making the range really large. This is of course useful to see the entire picture. So the finemerge is just in addition to this.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not make it so that if the dynamic range of the plot is too large it duplicates the plot and makes a "full range" version and a "zoomed in" version? Might be hard, but would make more sense from an ease of use and automation standpoint.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is hard to define 'too large'. Actually I think it would be best that we always have such 1.05-0.95 range canvas. So the too large would be defined as 'xmin<0.95 || xmax > 1.05'. And these days, it is almost always the case.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, then change the code to do that.

TString name = Form("ClosureVs%s_Overview_%s",getVariableTitleString(var).c_str(),alg.Data());
if(!flavor.IsNull()) name+="_"+flavor;

if (finemerge) name+="_Fine";
//
// Setup the frame, canvas, legend, and pave
//
Expand All @@ -548,10 +549,11 @@ void ClosureMaker::makeMergedCanvas() {
frame->GetXaxis()->SetLimits(XminCalo[0],Xmax[0]);
frame->GetXaxis()->SetMoreLogLabels();
frame->GetXaxis()->SetNoExponent();
pair<double,double> range = determineCanvasRange(frame->GetXaxis()->GetXmin(),frame->GetXaxis()->GetXmax());
//frame->GetYaxis()->SetRangeUser(0.95,1.05);
pair<double,double> range = determineCanvasRange(frame->GetXaxis()->GetXmin(),frame->GetXaxis()->GetXmax());
if (range.second - range.first > 0.2 && !finemerge) makeMergedCanvas(true);
if (finemerge) frame->GetYaxis()->SetRangeUser(0.95,1.05);
//frame->GetYaxis()->SetRangeUser(0.35,1.35);
frame->GetYaxis()->SetRangeUser(range.first,range.second);
else frame->GetYaxis()->SetRangeUser(range.first,range.second);
frame->GetXaxis()->SetTitle(getVariableAxisTitleString(var).c_str());
frame->GetYaxis()->SetTitle("Response");
canvases_legends.push_back(make_pair(tdrCanvas(name,frame,14,11,true),
Expand Down Expand Up @@ -635,5 +637,3 @@ void ClosureMaker::makeClosure(const VARIABLES::Variable ivar) {
loopOverAlgorithms();
closeFiles();
}


Loading