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

Update L2Creator.cc #21

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Update L2Creator.cc #21

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
168 changes: 157 additions & 11 deletions JetUtilities/src/L2Creator.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ L2Creator::L2Creator() {
histogramMetric = HistUtil::getHistogramMetricType(histMet);
delphes = false;
maxFitIter = 30;
ptclipfit=false;
}

//______________________________________________________________________________
Expand All @@ -46,6 +47,10 @@ L2Creator::L2Creator(CommandLine& cl) {
histMet = cl.getValue<string> ("histMet", "mu_h");
histogramMetric = HistUtil::getHistogramMetricType(histMet);

ptclip = cl.getValue<float> ("ptclip", 0.);
statTh = cl.getValue<int> ("statTh", 4);
ptclipfit = cl.getValue<bool> ("ptclipfit", false);

if (!cl.partialCheck()) return;
cl.print();
}
Expand Down Expand Up @@ -211,6 +216,14 @@ void L2Creator::loopOverEtaBins() {
TH1F* hrsp(0);
hl_rsp.begin_loop();

//print fit results for all eta bins in a txt file
TString txtFitResultsFilename = outputDir + "Fit_Results.txt";
ofstream FitResults(txtFitResultsFilename);
FitResults.setf(ios::right);

//edw gia draw prob,chi2 vs eta
double etabin[82], Chi2NDF[82], Prob[82], Prob_0p009[82];

while ((hrsp=hl_rsp.next_object(indices))) {

unsigned int ieta=indices[0];
Expand All @@ -232,7 +245,7 @@ void L2Creator::loopOverEtaBins() {
// only add points to the graphs if the current histo is not empty
// the current setting might be a little high
//
if (hrsp->GetEntries() > 4) {//hrsp->Integral()!=0) {
if (hrsp->GetEntries() > 30) {//hrsp->Integral()!=0) {//EDW 4

//TF1* frsp = (TF1*)hrsp->GetListOfFunctions()->Last();
//std::cout << "hrspName = " << hrsp->GetName() << ": frsp = " << frsp << std::endl;
Expand Down Expand Up @@ -277,19 +290,23 @@ void L2Creator::loopOverEtaBins() {
double eabsrsp = epeak;
double abscor = 0.0;
double eabscor = 0.0;

if (absrsp > 0) {
//cout << "EDW 1 " << "eta" << vabscor_eta.back()->GetName() << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl;
if (absrsp > 0 ) {//edw test
abscor =1.0/absrsp;
eabscor = abscor*abscor*epeak;


//cout << "EDW 2 " << "eta" << vabscor_eta.back()->GetName() << endl;
}
if ((abscor>0) && (absrsp>0) && (eabscor>1e-5) && (eabscor/abscor<0.5) && (eabsrsp>1e-4) && (eabsrsp/absrsp<0.5)) {
if ((abscor>0) && (absrsp>0) && (eabscor>1e-5) && (eabscor/abscor<0.5) && (eabsrsp>1e-4) && (eabsrsp/absrsp<0.5)) {
int n = vabsrsp_eta.back()->GetN();
vabsrsp_eta.back()->SetPoint (n,refpt, absrsp);
vabsrsp_eta.back()->SetPointError(n,erefpt,eabsrsp);
vabscor_eta.back()->SetPoint (n,jetpt, abscor);
vabscor_eta.back()->SetPointError(n,ejetpt,eabscor);
// cout << "EDW 3" << "eta" << vabscor_eta.back()->GetName() << "refpt" << refpt << "jetpt" << jetpt << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl;
}
else cout << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl;
//else cout << "EDW 4" << "absrsp " << absrsp << " and eabsrsp " << eabsrsp << " and abscor " << abscor << " and eabscor " << eabscor << endl;
}

//
Expand Down Expand Up @@ -376,6 +393,25 @@ void L2Creator::loopOverEtaBins() {
fabscor=new TF1("fit",fcn.Data(),xmin,xmax);
setOfflinePFParameters(gabscor, fabscor,xmin,xmax);

std::cout<<"ieta = "<<ieta<<std::endl;
// if(ieta==14 || ieta==67) std::cout<<"eta bin = "<<vabscor_eta.back()->GetName()<<std::endl;

/* if(ieta<=14 || ieta>=67) //EDW Mikko's fit function, fix parameters 3,4,5 to 0 for |eta|>=2.5
{
// fabscor->SetRange(12.,3500.); //for mixed fit range
fabscor->FixParameter(3,0.);
fabscor->FixParameter(4,0.);
fabscor->FixParameter(5,0.);
}
*/

if(alg.find("puppi")!=string::npos){
if(ieta<4 || ieta>77){ //for PUPPI
fabscor->FixParameter(3,0.);
fabscor->FixParameter(6,0.);
}
}

if(l2pffit.Contains("spline",TString::kIgnoreCase)) {
vabscor_eta_spline.back()->setPartialFunction(fabscor);
}
Expand Down Expand Up @@ -537,6 +573,54 @@ void L2Creator::loopOverEtaBins() {
}
}


//edw ptclipfit
if (ptclipfit)
{
if (xmin > 0.0001)
{
int nPar = fabscor->GetNpar();
int clipPar = nPar;

TString clip = TString::Format("((x<10)*([%d]))+((x>=10)*("+(TString)fabscor->GetTitle()+"))", clipPar);
TF1 * fabscornew = new TF1(fabscor->GetName(), clip, 0.001, 6500);
for (int ip=0; ip<nPar; ip++)
{
fabscornew->SetParameter(ip, fabscor->GetParameter(ip));
}

fabscornew->SetParameter(clipPar, fabscor->Eval(xmin));

fabscornew->SetChisquare(fabscor->GetChisquare());
fabscornew->SetNDF(fabscor->GetNDF());

fabscor = fabscornew;
gabscor->GetListOfFunctions()->Clear();
gabscor->GetListOfFunctions()->AddLast(fabscor);
}
}

//EDW print chi2 and prob for each fit in every eta bin
std::cout<<"Chi2/NDF = "<<fabscor->GetChisquare()/fabscor->GetNDF()<<std::endl;
std::cout<<"Prob = "<<fabscor->GetProb()<<std::endl;

FitResults<<ieta<<"\n";
FitResults<<"Eta bin : "<<vabscor_eta.back()->GetName()<<"\n";
FitResults<<"Chi2/NDF = "<<fabscor->GetChisquare()<<"/"<<fabscor->GetNDF()<<" = "<<fabscor->GetChisquare()/fabscor->GetNDF()<<"\n";
FitResults<<"Prob = "<<fabscor->GetProb()<<"\n\n\n";
FitResults<<"Number of points = "<<npoints<<"\n\n\n";


//===== for drawing graph Chi2/NDF vs Eta =====
etabin[ieta] = (hl_rsp.minimum(0,ieta)+hl_rsp.maximum(0,ieta))/2;
if(fabscor->GetChisquare()/fabscor->GetNDF()>1000) {Chi2NDF[ieta] =0; cout<<ieta<<endl;}
else Chi2NDF[ieta] = fabscor->GetChisquare()/fabscor->GetNDF();

Prob[ieta] = fabscor->GetProb();
if(fabscor->GetProb()<0.009) Prob_0p009[ieta] = fabscor->GetProb();
else Prob_0p009[ieta] = 0.01;


//
// format the graphs
//
Expand All @@ -553,6 +637,36 @@ void L2Creator::loopOverEtaBins() {
vabscor_eta_spline.push_back(nullptr);
}
}

FitResults.close();


//===== Draw graph Chi2/NDF vs Eta =====
TGraph *g_chi2ndf_vs_eta = new TGraph(82, etabin, Chi2NDF);
TCanvas *pad1 = new TCanvas("pad1", "",1);
g_chi2ndf_vs_eta->GetYaxis()->SetTitle("Chi2/ndf");
g_chi2ndf_vs_eta->GetXaxis()->SetTitle("eta");
g_chi2ndf_vs_eta->SetMarkerStyle(20);
g_chi2ndf_vs_eta->Draw("AP");
pad1->SaveAs(Form("%s/Chi2NDF_vs_eta.png",outputDir.Data()));

TGraph *g_prob_vs_eta = new TGraph(82, etabin, Prob);
TCanvas *pad2 = new TCanvas("pad2", "",1);
g_prob_vs_eta->GetYaxis()->SetTitle("Prob.");
g_prob_vs_eta->GetXaxis()->SetTitle("eta");
g_prob_vs_eta->SetMarkerStyle(20);
g_prob_vs_eta->Draw("AP");
pad2->SaveAs(Form("%s/Prob_vs_eta.png",outputDir.Data()));

TGraph *g_prob_0p009_vs_eta = new TGraph(82, etabin, Prob_0p009);
TCanvas *pad3 = new TCanvas("pad3", "",1);
g_prob_0p009_vs_eta->GetYaxis()->SetTitle("Prob");
g_prob_0p009_vs_eta->GetXaxis()->SetTitle("eta");
g_prob_0p009_vs_eta->SetMarkerStyle(20);
g_prob_0p009_vs_eta->Draw("AP");
pad3->SaveAs(Form("%s/prob_0p009_vs_eta.png",outputDir.Data()));


}

//______________________________________________________________________________
Expand Down Expand Up @@ -804,7 +918,7 @@ void L2Creator::makeCanvas(string makeCanvasVariable) {
spline_func->SetRange(bounds.first,bounds.second);
spline_func->SetLineColor(igraph%nperpad+1);
spline_func->SetLineStyle(kDotted);
spline_func->Draw("same");
// spline_func->Draw("same");
}
}

Expand Down Expand Up @@ -874,7 +988,8 @@ TString L2Creator::getOfflinePFFunction() {
return "[0]+([1]/(pow(log10(x),2)+[2]))+([3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5])))";
}
else if(l2pffit.EqualTo("standard+Gaussian",TString::kIgnoreCase)) {
return "[0]+([1]/(pow(log10(x),2)+[2]))+([3]*exp(-([4]*((log10(x)-[5])*(log10(x)-[5])))))+([6]*exp(-([7]*((log10(x)-[8])*(log10(x)-[8])))))";
return "[0]+([1]/(pow(log10(x),2)+[2]))+([3]*exp(-([4]*((log10(x)-[5])*(log10(x)-[5])))))+([6]*exp(-([7]*((log10(x)-[8])*(log10(x)-[8])))))";
// return "[0]+[1]*pow(x,[2])+[3]*exp(-[4]*max(0.,(log10(x)-[5]))*(log10(x)-[5]))"; //Modification of 0,1,2 with powerlaw first term + single-sided log-gaus by Mikko
}
else if(l2pffit.EqualTo("LogNormal+Gaussian+fixed",TString::kIgnoreCase)) {
return "([0]+TMath::LogNormal(TMath::Log10(x),[1],[2],[3]))+([4]+[5]/(pow(log10(x),2)+[6])+[7]*exp(-[8]*(log10(x)-[9])*(log10(x)-[9])))";
Expand Down Expand Up @@ -939,7 +1054,7 @@ void L2Creator::setOfflinePFParameters(TGraphErrors* gabscor, TF1* fabscor, doub
fabscor->SetParLimits(4,0,100);
}
else if(l2pffit.EqualTo("standard+Gaussian",TString::kIgnoreCase)) {
fabscor->SetParameter(0,-0.0221278);
/* fabscor->SetParameter(0,-0.0221278);
fabscor->SetParameter(1,119.265);
fabscor->SetParameter(2,100);
fabscor->SetParameter(3,-0.0679365);
Expand All @@ -954,6 +1069,34 @@ void L2Creator::setOfflinePFParameters(TGraphErrors* gabscor, TF1* fabscor, doub
fabscor->SetParLimits(4,0,500);
fabscor->SetParLimits(0,-2,25);
fabscor->SetParLimits(1,0,250);
*/
fabscor->SetParameter(0,0.0221278); //v2
fabscor->SetParameter(1,14.265);
fabscor->SetParameter(2,10);
fabscor->SetParameter(3,-0.0679365);
fabscor->SetParameter(4,2.82597);
fabscor->SetParameter(5,1.8277);
fabscor->SetParameter(6,-0.0679365);
fabscor->SetParameter(7,3.82597);
fabscor->SetParameter(8,1.8277);
fabscor->SetParLimits(6,-20,10);
fabscor->SetParLimits(7,0,100);
fabscor->SetParLimits(3,-15,15);
fabscor->SetParLimits(4,0,5);
fabscor->SetParLimits(5,0,50);
fabscor->SetParLimits(0,-2,50);
fabscor->SetParLimits(2,0,200);
fabscor->SetParLimits(1,0,250);


/* fabscor->SetRange(8.,3500.);
fabscor->SetParameter(0,1.09892); //Modification of 0,1,2 with powerlaw first term + single-sided log-gaus by Mikko
fabscor->SetParameter(1,-0.46266);
fabscor->SetParameter(2,-0.152203);
fabscor->SetParameter(3,0.0951505);
fabscor->SetParameter(4,1.16479);
fabscor->SetParameter(5,1.53048);
*/
}
else if(l2pffit.EqualTo("LogNormal+Gaussian+fixed",TString::kIgnoreCase)) {
fabscor->SetRange(3,2000);
Expand Down Expand Up @@ -1313,6 +1456,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() {

//For eta-dependent spline clipping
int pt_limit = 70;
float pt_clip = ptclip;

unsigned int vector_size = 0;
vector_size = vabscor_eta.size();
Expand Down Expand Up @@ -1345,6 +1489,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() {

bool abovePtLimit = false;
bool lastLine = false;
bool firstline = true;

for(int isection=0; isection<spline->getNSections(); isection++) {
if(lastLine) continue;
Expand All @@ -1364,7 +1509,7 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() {
//if(isection==spline->getNSections()-1) {
// bounds.second = 6500;
//}

if(bounds.second < pt_clip) continue;
if(isection==spline->getNSections()-1) lastLine = true;
if(bounds.second >= pt_limit) {
abovePtLimit = true;
Expand All @@ -1373,16 +1518,17 @@ void L2Creator::writeTextFileForCurrentAlgorithm_spline() {

//For expediency of Summer16_25nsV5_MC do eta-dependent clipping
fout<<setw(8) <<etamin<<setw(8)<<etamax
<<setw(10)<<setprecision(6)<<(isection ? bounds.first : 0.001)
<<setw(10)<<setprecision(6)<<(firstline ? 0.001 : bounds.first)
<<setw(10)<<setprecision(6)<<(lastLine ? 6500 : bounds.second)
<<setw(6)<<(int)(spline->getNpar()+2) //Number of parameters + 2
<<setw(12)<<setprecision(8)<<bounds.first
<<setw(12)<<setprecision(8)<<(firstline ? pt_clip : bounds.first)
<<setw(12)<<setprecision(8)<<(abovePtLimit ? pt_limit : bounds.second);
TF1* spline_func = spline->setParameters(isection);
for(int p=0; p<spline->getNpar(); p++) {
fout<<setw(17)<<setprecision(10)<<spline_func->GetParameter(p);
}
fout<<endl;
firstline = false;
}
}
}
Expand Down