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.

Sign up for GitHub

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

Jump to bottom

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

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

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
4b79df2
clip at 8 for low pt
sifuluo Feb 2, 2018
d960cb0
Adding functions of clipping pt in lower pt and statistical threshold…
sifuluo Feb 5, 2018
691f9be
Fit included in HistogramUtility
sifuluo Feb 11, 2018
d30f828
Update run_JRA_cfg.py
sifuluo Feb 13, 2018
4e7ddc3
reversed run_JRA_cfg.py to default
sifuluo Feb 13, 2018
a388b43
Merge branch 'master' of github.com:sifuluo/JetMETAnalysis
sifuluo Feb 13, 2018
ca0d9b0
instance with nentries <= 4 will not be shown in closure
sifuluo Mar 5, 2018
b7b393d
updated
sifuluo Jul 2, 2018
d7e8847
minor update ready to merge.
sifuluo Jul 5, 2018
eee7200
compareJEC now works with commandline configurations and can generate…
sifuluo Aug 8, 2018
b408bdf
compareJEC now works with commandline configurations and can generate…
sifuluo Aug 8, 2018
eab6e92
In L2Creator, ptclip is kept while ptclips and ptclipcones is added t…
sifuluo Aug 8, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
187 changes: 88 additions & 99 deletions JetAnalyzers/bin/compareJEC.C
View file
Open in desktop

Large diffs are not rendered by default.

208 changes: 104 additions & 104 deletions JetAnalyzers/bin/jet_synchplot_x.cc
View file
Open in desktop

Large diffs are not rendered by default.

22 changes: 11 additions & 11 deletions JetAnalyzers/bin/jet_synchtest_x.cc
View file
Open in desktop
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 @@ -796,7 +797,6 @@ void MatchEventsAndJets::LoopOverEvents(bool verbose, bool reduceHistograms, str

if (iftest && nevs >= maxEvts) return;

//if (nevs%10000==0) cout << "\t"<<nevs << endl;
loadbar2(nevs+1,nentries,50,"\t\t");

// if this entry does not exist on the second ntuple just skip this event
Expand Down Expand Up @@ -830,7 +830,7 @@ void MatchEventsAndJets::LoopOverEvents(bool verbose, bool reduceHistograms, str
else {
jetMapIndex++;
ReadJetMap(jetMapIndex,readJetMap);
}
}

if(FillHistograms(reduceHistograms)) nevs++;

Expand Down Expand Up @@ -908,7 +908,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 +1186,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 +1356,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 +1372,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 +1389,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
View file
Open in desktop
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
View file
Open in desktop
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>
80 changes: 40 additions & 40 deletions JetUtilities/interface/ClosureMaker.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,48 +62,48 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
class ClosureMaker {
public:
ClosureMaker();
ClosureMaker(CommandLine& cl);

//
// member functions
//
void openInputFile();
void getHistograms(TDirectoryFile* idir);
void openOutputFile();
void closeFiles();
void makeLines();
void loopOverDirectories();
void loopOverAlgorithms();
void resetForNextAlgorithm();
void loopOverHistograms();
void etaClosureMergedPtBins(TDirectoryFile* idir);
void loopOverBins(TH2F* hvar, unsigned int iVarBin);
void adjust_fitrange(TH1* h,double& min,double& max);
void checkResponse();
pair<double,double> determineCanvasRange(double xmin, double xmax);
void makeCanvases();
void makeMergedCanvas();
void writeToFile();
void makeClosure(const VARIABLES::Variable ivar = VARIABLES::refpt);
ClosureMaker();
ClosureMaker(CommandLine& cl);
//
// member functions
//
void openInputFile();
void getHistograms(TDirectoryFile* idir);
void openOutputFile();
void closeFiles();
void makeLines();
void loopOverDirectories();
void loopOverAlgorithms();
void resetForNextAlgorithm();
void loopOverHistograms();
void etaClosureMergedPtBins(TDirectoryFile* idir);
void loopOverBins(TH2F* hvar, unsigned int iVarBin);
void adjust_fitrange(TH1* h,double& min,double& max);
void checkResponse();
pair<double,double> determineCanvasRange(double xmin, double xmax);
void makeCanvases();
void makeMergedCanvas(bool finemerge);
void writeToFile();
void makeClosure(const VARIABLES::Variable ivar = VARIABLES::refpt);

private:
bool objects_loaded, draw_guidelines;
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<TH1F*> hClosure;
TF1 *line, *linePlus, *lineMinus;
vector<pair<TCanvas*,TLegend*> > canvases_legends;
vector<TPaveText*> pave;
VARIABLES::Variable var;
TDirectoryFile *odir;
HistUtil::HistogramMetric histogramMetric;
bool objects_loaded, draw_guidelines;
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<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
View file
Open in desktop
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
5 changes: 4 additions & 1 deletion JetUtilities/interface/L2Creator.hh
View file
Open in desktop
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, statThreshold;
HistUtil::HistogramMetric histogramMetric;
TFile* ofile;
TFile* ifile;
Expand All @@ -116,6 +116,9 @@ private:
vector<TGraphErrors*> vabscor_eta;
vector<TGraph*> vrelcor_eta;
vector<PiecewiseSpline*> vabscor_eta_spline;
vector<int> ptclipcones;
vector<float> ptclips;
float ptclip;
};

#endif
Loading