本文整理汇总了C++中TGraph::Eval方法的典型用法代码示例。如果您正苦于以下问题:C++ TGraph::Eval方法的具体用法?C++ TGraph::Eval怎么用?C++ TGraph::Eval使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TGraph
的用法示例。
在下文中一共展示了TGraph::Eval方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: getValsFromLikelihood
vector<double> getValsFromLikelihood(TGraph *graph){
TGraph *grRot = new TGraph();
TGraph *grRotUpp = new TGraph();
TGraph *grRotLow = new TGraph();
int pLow=0,pHigh=0;
pair<double,double> graphMin = getGraphMin(graph);
for (int p=0; p<graph->GetN(); p++){
double x,y;
graph->GetPoint(p,x,y);
grRot->SetPoint(p,y,x);
if (x<=graphMin.first) {
grRotLow->SetPoint(pLow,y,x);
pLow++;
}
if (x>=graphMin.first) {
grRotUpp->SetPoint(pHigh,y,x);
pHigh++;
}
}
// return best fit and +/- 1/2 sigma errors
vector<double> result;
result.push_back(grRot->Eval(0.));
result.push_back(grRotLow->Eval(1.0));
result.push_back(grRotUpp->Eval(1.0));
result.push_back(grRotLow->Eval(4.0));
result.push_back(grRotUpp->Eval(4.0));
return result;
}
示例2: points
void points(TString filename) {
TString cmssw;
// 167
cmssw = "$1.6.7$";
TFile *f = TFile::Open(filename);
std::vector< TString > taggers;
taggers.push_back( "gTC2_udsg" );
taggers.push_back( "gTC3_udsg" );
taggers.push_back( "gTP_udsg" );
taggers.push_back( "gJBP_udsg" );
taggers.push_back( "gSSV_udsg" );
taggers.push_back( "gCSV_udsg" );
std::vector< TString > discriminators;
discriminators.push_back( "discTC2_udsg" );
discriminators.push_back( "discTC3_udsg" );
discriminators.push_back( "discTP_udsg" );
discriminators.push_back( "discJBP_udsg" );
discriminators.push_back( "discSSV_udsg" );
discriminators.push_back( "discCSV_udsg" );
//TCanvas *cv_TC = new TCanvas("cv_TC","cv_TC",700,700);
//TCanvas *cv_TP = new TCanvas("cv_TP","cv_TP",700,700);
std::cout << "Tagger & Point & Discriminator & light mistag & b-efficiency \\\\ \\hline" << std::endl;
for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
TString tag = taggers[itagger];
TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("Histograms/MCtruth/"+tag);
//if (taggers == "gTC2_udsg" || taggers =
TGraph *dgraph = (TGraph*) gDirectory->Get("Histograms/MCtruth/"+discriminators[itagger]);
dgraph->Sort();
TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
g->Sort();
//cv[itagger] = new TCanvas("cv","cv",600,600);
//g->Draw("ACP");
TString se = " & ";
std::cout << tag << se << "Loose" << se << std::setprecision(3) << dgraph->Eval(0.1) << se << "0.1" << se << std::setprecision(2) << g->Eval(0.1) << "\\\\" << std::endl;
std::cout << tag << se << "Medium" << se << std::setprecision(3) << dgraph->Eval(0.01) << se << "0.01" << se << std::setprecision(2) << g->Eval(0.01) << "\\\\" << std::endl;
std::cout << tag << se << "Tight" << se << std::setprecision(3) << dgraph->Eval(0.001) << se << "0.001" << se << std::setprecision(2) << g->Eval(0.001) << "\\\\ \\hline" << std::endl;
}
}
示例3: tit
/**
* Create ratios to other data
*
* @param ib Bin number
* @param res Result
* @param alice ALICE result if any
* @param cms CMS result if any
* @param all Stack to add ratio to
*/
void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
{
if (!all || !res || !(alice || cms)) return;
Int_t off = 5*ib;
TGraph* gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
TGraph** pg = gs;
while (*pg) {
TGraph* g = *pg;
const char* n = (g == alice ? "ALICE" : "CMS");
TH1* r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
TString tit(r->GetTitle());
tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
r->SetTitle(tit);
r->SetMarkerColor(g->GetMarkerColor());
r->SetLineColor(g->GetLineColor());
TObject* tst = r->FindObject("legend");
if (tst) r->GetListOfFunctions()->Remove(tst);
for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
Double_t c = r->GetBinContent(i);
Double_t e = r->GetBinError(i);
Double_t o = g->Eval(r->GetBinCenter(i));
if (o < 1e-12) {
r->SetBinContent(i, 0);
r->SetBinError(i, 0);
continue;
}
r->SetBinContent(i, (c - o) / o + off);
r->SetBinError(i, e / o);
}
all->Add(r);
pg++;
}
TLegend* leg = StackLegend(all);
if (!leg) return;
TString txt = res->GetTitle();
txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
if (ib == 0) txt.Append(" "); // (#times1)");
// else if (ib == 1) txt.Append(" (#times10)");
else txt.Append(Form(" (+%d)", off));
TObject* dummy = 0;
TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
e->SetMarkerStyle(res->GetMarkerStyle());
e->SetMarkerSize(res->GetMarkerSize());
e->SetMarkerColor(kBlack);
e->SetFillColor(0);
e->SetFillStyle(0);
e->SetLineColor(kBlack);
}
示例4: disceff
void disceff(TString filename) {
gROOT->SetStyle("Plain");
TString cmssw;
// 167
cmssw = "$3.1.0_pre9$";
TFile *f = TFile::Open(filename);
std::vector< TString > taggers;
taggers.push_back( "TC2" );
taggers.push_back( "TC3" );
taggers.push_back( "TP" );
taggers.push_back( "BTP" );
taggers.push_back( "SSV" );
taggers.push_back( "CSV" );
taggers.push_back( "MSV" );
taggers.push_back( "SMT" );
taggers.push_back( "SETbyIP3d" );
taggers.push_back( "SETbyPt" );
taggers.push_back( "SMTbyIP3d" );
taggers.push_back( "SMTbyPt" );
std::vector< TString > discriminators;
for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
discriminators.push_back( "disc"+taggers[itagger]+"_udsg" );
}
// discriminators.push_back( "discTC3_udsg" );
// discriminators.push_back( "discTP_udsg" );
const int dim=taggers.size();
TCanvas *cv[dim];
TMultiGraph* mg[dim];
for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
TString tag = "g"+taggers[itagger]+"_udsg";
TString tagb = "g"+taggers[itagger]+"_b";
TString tagc = "g"+taggers[itagger]+"_c";
cv[itagger] = new TCanvas("cv_"+taggers[itagger],"cv_"+taggers[itagger],700,700);
TLegend *legend0 = new TLegend(0.68,0.70,0.88,0.90);
TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tag);
TGraphErrors *bgraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tagb);
TGraphErrors *cgraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tagc);
TGraph *dgraph = (TGraph*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+discriminators[itagger]);
TGraph *bvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),bgraph->GetY());
TGraph *cvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),cgraph->GetY());
TGraph *lvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),agraph->GetY());
TGraph *udsgvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),dgraph->GetX());
dgraph->Sort();
// udsgvsdgraph->Sort();
// udsgvsdgraph->SetLineColor(1);
// legend0 -> AddEntry(udsgvsdgraph,"control","l");
lvsdgraph->Sort();
lvsdgraph->SetLineColor(2);
legend0 -> AddEntry(lvsdgraph,tag,"l");
cvsdgraph->Sort();
cvsdgraph->SetLineColor(3);
legend0 -> AddEntry(cvsdgraph,tagc,"l");
bvsdgraph->Sort();
bvsdgraph->SetLineColor(4);
legend0 -> AddEntry(bvsdgraph,tagb,"l");
mg[itagger]= new TMultiGraph();
// mg[itagger]->Add(udsgvsdgraph);
mg[itagger]->Add(lvsdgraph);
mg[itagger]->Add(cvsdgraph);
mg[itagger]->Add(bvsdgraph);
// mg[itagger]->Add(dgraph);
cv[itagger]->cd(1);
mg[itagger]->Draw("ALP");
mg[itagger]->GetYaxis()->SetTitle("eff");
mg[itagger]->GetXaxis()->SetTitle("discriminant");
legend0 -> Draw();
cv[itagger]->Update();
cv[itagger]->cd(0);
cv[itagger]-> Print ("BTagPATeff_vs_disc"+taggers[itagger]+".eps");
cv[itagger]-> Print ("BTagPATeff_vs_disc"+taggers[itagger]+".ps");
TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
g->Sort();
g->SetLineColor(itagger+1);
std::cout << " Tagger: " << tag << std::endl;
std::cout << " Loose(10%): " << " cut > " << std::setprecision(4) << dgraph->Eval(0.1) << " b-eff = " << g->Eval(0.1) << std::endl;
std::cout << " Medium(1%): " << " cut > " << std::setprecision(4) << dgraph->Eval(0.01) << " b-eff = " << g->Eval(0.01) << std::endl;
std::cout << " Tight(0.1%) = " << " cut > " << std::setprecision(4) << dgraph->Eval(0.001) << " b-eff = " << g->Eval(0.001) << std::endl;
}//end for
}
示例5: DMplot
//.........这里部分代码省略.........
h_scalar_min->GetXaxis()->SetTitleOffset(1.03);
h_scalar_min->GetXaxis()->SetTitle("DM mass [GeV]");
h_scalar_min->GetYaxis()->SetTitle("DM-nucleon cross section [cm^{2}]");
h_scalar_lat->SetLineColor(4);
h_scalar_lat->SetLineStyle(1);
h_scalar_lat->SetLineWidth(3);
h_scalar_max->SetLineColor(4);
h_scalar_max->SetLineStyle(2);
h_scalar_max->SetLineWidth(3);
h_fermion_min->SetLineColor(kRed);
h_fermion_min->SetLineStyle(2);
h_fermion_min->SetLineWidth(3);
h_fermion_lat->SetLineColor(kRed);
h_fermion_lat->SetLineStyle(1);
h_fermion_lat->SetLineWidth(3);
h_fermion_max->SetLineColor(kRed);
h_fermion_max->SetLineStyle(2);
h_fermion_max->SetLineWidth(3);
//h_scalar_lat->SetFillStyle(3005);
//h_scalar_lat->SetLineWidth(-402);
h_scalar_min->Draw("AL");
h_scalar_lat->Draw("L");
h_scalar_max->Draw("L");
h_fermion_min->Draw("L");
h_fermion_lat->Draw("L");
h_fermion_max->Draw("L");
z12016->Draw("L");
//z12015->Draw("L");
//CRESST2->Draw("C");
SCDMS->Draw("C");
PANDAX->Draw("C");
//leg2->Draw();
TLatex *lat = new TLatex(); lat->SetTextSize(0.025);
lat->SetTextFont(42);
lat->SetTextColor(kGreen+2);
lat->SetTextAngle(15);
//lat->DrawLatex(130,z1->Eval(130)*1.5,"LUX #it{Phys. Rev. Lett.} #bf{116} (2016)");
//lat->DrawLatex(130,z12015->Eval(130)*1.5,"LUX (2015)");
lat->DrawLatex(130,z12016->Eval(130)*0.5,"LUX (2013+2014-16)");
lat->SetTextColor(SCDMS->GetLineColor());
lat->SetTextAngle(344);
//lat->SetFillColor(kWhite);
lat->DrawLatex(5,SCDMS->Eval(5)*1.5,"CDMSlite (2015)");
lat->SetTextColor(kBlue);
lat->SetTextAngle(330);
lat->DrawLatex(13,h_scalar_lat->Eval(13)*1.5,"Scalar DM");
lat->SetTextColor(kRed);
lat->SetTextAngle(4);
lat->DrawLatex(5,h_fermion_lat->Eval(5)*1.5,"Fermion DM");
lat->SetTextColor(kMagenta+2);
lat->SetTextAngle(15);
lat->DrawLatex(130,PANDAX->Eval(130)*1.25,"PandaX-II (2016)");
TLatex * tex = new TLatex();
tex->SetNDC();
tex->SetTextFont(42);
tex->SetLineWidth(2);
tex->SetTextSize(0.04);
tex->SetTextAlign(31);
tex->DrawLatex(0.93,0.78,"4.9 fb^{-1} (7 TeV) + 19.7 fb^{-1} (8 TeV)");
tex->DrawLatex(0.93,0.74,"+ 2.3 fb^{-1} (13 TeV)");
//tex->SetTextAlign(11);
tex->SetTextFont(42);
tex->SetTextSize(0.06);
TLegend *legBOX;
if (isPrelim) {
legBOX = new TLegend(0.16,0.83,0.38,0.89);
legBOX->SetFillColor(kWhite);
legBOX->SetLineWidth(0);
legBOX->Draw();
tex->DrawLatex(0.17, 0.84, "#bf{CMS} #it{Preliminary}");
} else {
legBOX = new TLegend(0.16,0.83,0.28,0.89);
legBOX->SetFillColor(kWhite);
legBOX->SetLineWidth(0);
//legBOX->Draw();
tex->DrawLatex(0.93, 0.84, "#bf{CMS}");
}
tex->SetTextSize(0.045);
tex->DrawLatex(0.93,0.68,Form("B(H #rightarrow inv.) < %.2f",BRinv));
tex->DrawLatex(0.93,0.5,"90% CL limits");
canv->SetTicky(1);
canv->SetTickx(1);
canv->SaveAs("limitsDM.pdf");
}
示例6: mass4Chan
//.........这里部分代码省略.........
graph->SetPoint(182,94.90000153,1.380649686);
graph->SetPoint(183,94.94000244,1.434825897);
graph->SetPoint(184,94.98000336,1.490056992);
graph->SetPoint(185,95.01999664,1.546342969);
graph->SetPoint(186,95.05999756,1.603683949);
graph->SetPoint(187,95.09999847,1.66207993);
graph->SetPoint(188,95.13999939,1.721530676);
graph->SetPoint(189,95.18000031,1.782036185);
graph->SetPoint(190,95.22000122,1.843595982);
graph->SetPoint(191,95.22000122,1.843595982);
graph->SetPoint(192,95.26000214,1.906209826);
graph->SetPoint(193,95.30000305,1.969877124);
graph->SetPoint(194,95.33999634,2.034597158);
graph->SetPoint(195,95.37999725,2.100369453);
graph->SetPoint(196,95.41999817,2.167192936);
graph->SetPoint(197,95.45999908,2.235066652);
graph->SetPoint(198,95.5,2.30398941);
graph->SetPoint(199,95.54000092,2.37395978);
graph->SetPoint(200,95.58000183,2.44497633);
graph->SetPoint(201,95.62000275,2.517037392);
graph->SetPoint(202,95.66000366,2.590141296);
graph->SetPoint(203,95.69999695,2.66428566);
graph->SetPoint(204,95.73999786,2.739468575);
graph->SetPoint(205,95.77999878,2.815687418);
graph->SetPoint(206,95.81999969,2.892939568);
graph->SetPoint(207,95.86000061,2.971222401);
graph->SetPoint(208,95.90000153,3.050532341);
graph->SetPoint(209,95.94000244,3.130866289);
graph->SetPoint(210,95.98000336,3.212220907);
graph->Sort();
cout << "4e: " << graph->Eval(91.1876) << endl;
//2e2mu
TGraph *graph1 = new TGraph(211);
graph1->SetName("Graph1");
graph1->SetTitle("Graph1");
graph1->SetFillColor(1);
graph1->SetLineWidth(3);
graph1->SetLineColor(kBlue);
graph1->SetMarkerStyle(20);
graph1->SetPoint(0,88.01999664,8.795412064);
graph1->SetPoint(1,88.05999756,8.61370182);
graph1->SetPoint(2,88.09999847,8.43285656);
graph1->SetPoint(3,88.13999939,8.252916336);
graph1->SetPoint(4,88.18000031,8.073918343);
graph1->SetPoint(5,88.22000122,7.895903111);
graph1->SetPoint(6,88.26000214,7.718908787);
graph1->SetPoint(7,88.30000305,7.542974472);
graph1->SetPoint(8,88.33999634,7.36813879);
graph1->SetPoint(9,88.37999725,7.194441795);
graph1->SetPoint(10,88.41999817,7.021921158);
graph1->SetPoint(11,88.45999908,6.850616932);
graph1->SetPoint(12,88.5,6.680567265);
graph1->SetPoint(13,88.54000092,6.51181221);
graph1->SetPoint(14,88.58000183,6.344389915);
graph1->SetPoint(15,88.62000275,6.178339481);
graph1->SetPoint(16,88.66000366,6.013700008);
graph1->SetPoint(17,88.69999695,5.850510597);
graph1->SetPoint(18,88.73999786,5.688809872);
graph1->SetPoint(19,88.77999878,5.528635979);
graph1->SetPoint(20,88.81999969,5.370028496);
graph1->SetPoint(21,88.81999969,5.370028496);
示例7: plotLimit
void plotLimit(string outputDir="./", TString inputs="", TString inputs_blinded="", TString inputXSec="", bool strengthLimit=true, bool blind=false, double energy=7, double luminosity=5.035, TString legendName="ee and #mu#mu channels")
{
setTDRStyle();
gStyle->SetPadTopMargin (0.05);
gStyle->SetPadBottomMargin(0.12);
gStyle->SetPadRightMargin (0.16);
gStyle->SetPadLeftMargin (0.14);
gStyle->SetTitleSize(0.04, "XYZ");
gStyle->SetTitleXOffset(1.1);
gStyle->SetTitleYOffset(1.45);
gStyle->SetPalette(1);
gStyle->SetNdivisions(505);
//get the limits from the tree
TFile* file = TFile::Open(inputs);
printf("Looping on %s\n",inputs.Data());
if(!file) return;
if(file->IsZombie()) return;
TFile* file_blinded = TFile::Open(inputs_blinded);
printf("Looping on %s\n",inputs_blinded.Data());
if(!file_blinded) return;
if(file_blinded->IsZombie()) return;
TTree* tree_blinded = (TTree*)file_blinded->Get("limit");
tree_blinded->GetBranch("mh" )->SetAddress(&Tmh );
tree_blinded->GetBranch("limit" )->SetAddress(&Tlimit );
tree_blinded->GetBranch("limitErr" )->SetAddress(&TlimitErr);
tree_blinded->GetBranch("quantileExpected")->SetAddress(&TquantExp);
TGraph* ExpLimitm2 = getLimitGraph(tree_blinded,0.025);
TGraph* ExpLimitm1 = getLimitGraph(tree_blinded,0.160);
TGraph* ExpLimit = getLimitGraph(tree_blinded,0.500);
TGraph* ExpLimitp1 = getLimitGraph(tree_blinded,0.840);
TGraph* ExpLimitp2 = getLimitGraph(tree_blinded,0.975);
file_blinded->Close();
TTree* tree = (TTree*)file->Get("limit");
tree->GetBranch("mh" )->SetAddress(&Tmh );
tree->GetBranch("limit" )->SetAddress(&Tlimit );
tree->GetBranch("limitErr" )->SetAddress(&TlimitErr);
tree->GetBranch("quantileExpected")->SetAddress(&TquantExp);
TGraph* ObsLimit = getLimitGraph(tree,-1 );
file->Close();
FILE* pFileSStrenght = fopen((outputDir+"SignalStrenght").c_str(),"w");
std::cout << "Printing Signal Strenght" << std::endl;
for(int i=0;i<ExpLimit->GetN();i++){
double M = ExpLimit->GetX()[i];
std::cout << "Mass: " << M << "; ExpLimit: " << ExpLimit->Eval(M) << std::endl;
printf("$%8.6E$ & $%8.6E$ & $[%8.6E,%8.6E]$ & $[%8.6E,%8.6E]$ \\\\\\hline\n",M, ExpLimit->Eval(M), ExpLimitm1->Eval(M), ExpLimitp1->Eval(M), ExpLimitm2->Eval(M), ExpLimitp2->Eval(M));
fprintf(pFileSStrenght, "$%8.6E$ & $%8.6E$ & $[%8.6E,%8.6E]$ & $[%8.6E,%8.6E]$ & $%8.6E$ \\\\\\hline\n",M, ExpLimit->Eval(M), ExpLimitm1->Eval(M), ExpLimitp1->Eval(M), ExpLimitm2->Eval(M), ExpLimitp2->Eval(M), ObsLimit->Eval(M));
if(int(ExpLimit->GetX()[i])%50!=0)continue; //printf("%f ",ObsLimit->Eval(M));
}printf("\n");
fclose(pFileSStrenght);
//get the pValue
inputs = inputs.ReplaceAll("/LimitTree", "/PValueTree");
file = TFile::Open(inputs);
printf("Looping on %s\n",inputs.Data());
if(!file) return;
if(file->IsZombie()) return;
tree = (TTree*)file->Get("limit");
tree->GetBranch("limit" )->SetAddress(&Tlimit );
TGraph* pValue = getLimitGraph(tree,-1);
file->Close();
//make TH Cross-sections
string suffix = outputDir;
TGraph* THXSec = Hxswg::utils::getXSec(outputDir);
scaleGraph(THXSec, 1000); //convert cross-section to fb
double cprime=1.0; double brnew=0.0;
double XSecScaleFactor = 1.0;
if(suffix.find("_cp")!=string::npos){
sscanf(suffix.c_str()+suffix.find("_cp"), "_cp%lf_brn%lf", &cprime, &brnew);
XSecScaleFactor = pow(cprime,2) * (1-brnew);
}
//XSecScaleFactor = 0.001; //pb to fb
scaleGraph(THXSec, XSecScaleFactor);
string prod = "pp_SM";
if(outputDir.find("ggH")!=std::string::npos)prod="gg";
if(outputDir.find("qqH")!=std::string::npos)prod="qq";
if(outputDir.find("ppH")!=std::string::npos)prod="pp";
strengthLimit = false;
if(prod=="pp_SM")strengthLimit=true;
//TGraph *XSecMELA = Hxswg::utils::getXSecMELA(cprime);
//Hxswg::utils::multiplyGraph( ObsLimit, XSecMELA);
//Hxswg::utils::multiplyGraph( ExpLimitm2, XSecMELA);
//Hxswg::utils::multiplyGraph( ExpLimitm1, XSecMELA);
//Hxswg::utils::multiplyGraph( ExpLimit, XSecMELA);
//Hxswg::utils::multiplyGraph( ExpLimitp1, XSecMELA);
//.........这里部分代码省略.........
示例8: AnalysisSparse
//.........这里部分代码省略.........
gr_true_eff->GetErrorX(i), gr_true_eff->GetErrorY(i));
if (!save_output) continue;
gSystem->mkdir(dir_prefix.Data());
tout = Form("%f %f %f %f", oux, ouy, ouxe, ouye);
if (i == 0)
tout = Form("Printf(\"%s\"); > %s/effi_%s", tout.Data(),
dir_prefix.Data(), grapht.Data());
else
tout = Form("Printf(\"%s\"); >> %s/effi_%s", tout.Data(),
dir_prefix.Data(), grapht.Data());
// Printf(":::::: %s", tout.Data());
gROOT->ProcessLine(tout.Data());
}
// ------------------
c = new TCanvas("cfinal", "mc_effi", 1200, 450);
c->Divide(2, 1, 0.001, 0.001); c->Modified(); c->Draw();
c->cd(1);
gr_true->SetMinimum(0);
gr_true->SetTitle(Form("phi (true & raw), %s", gtitle.Data()));
gr_true->SetMarkerColor(kGreen+1);
gr_true->Draw("AP");
gr->SetMarkerColor(kRed+1);
gr->Draw("P");
c->cd(2)->SetGrid();
gr_true_eff->SetMinimum(0);
gr_true_eff->SetTitle(Form("efficiency, %s", grapht.Data()));
gr_true_eff->SetMarkerColor(kGreen+1);
gr_true_eff->Draw("AP");
gr_e->SetMarkerColor(kRed+1);
gr_e->Draw("P");
// c->SaveAs(Form("%s_%s.eps", blabla.Data(), grapht.Data()));
return;
}
// TGraph *geff = new TGraph(Form("effi_raw_%s", grapht.Data()));
// TGraph *geff = new TGraph(Form("effi_true_%s", grapht.Data()));
// TGraph *geff = new TGraph("effi_true_Phi2010_qualityonly");
TGraph *geff = new TGraph("effi_true_PhiNsigma_qualityonly");
if (geff->IsZombie()) return;
geff->SetMarkerStyle(22);
geff->SetMarkerColor(kBlack);
geff->GetXaxis()->SetTitle("p_{t}, GeV/c");
geff->SetTitle(Form("efficiency, %s", grapht.Data()));
c = new TCanvas();
c->SetGrid();
geff->Draw("AP");
Double_t tpcsigma = 9999.9;
if (ilist == 1) tpcsigma = 1.0;
if (ilist == 2) tpcsigma = 1.5;
if (ilist == 3) tpcsigma = 2.0;
if (ilist == 4) tpcsigma = 2.5;
if (ilist == 5) tpcsigma = 3.0;
Double_t sss = TMath::Erf(tpcsigma/TMath::Sqrt(2.0));
if (noSigma) sss = 1.0;
Printf("sigma = %10f", sss);
// for (Int_t i = 0; i < count; i++)
// geff->GetY()[i] = (sss*sss)/(geff->GetY()[i]);
// geff->SetMaximum(1.0);
// geff->Draw("AP");
for (Int_t i = 0; i < count; i++) {
Double_t deno = geff->Eval(grx[i])*sss*sss;
if (deno < 0.00001) deno = 1;
gry_fix[i] = gry[i]/deno;
}
new TCanvas;
TGraph *gr_fix = new TGraph(count, grx, gry_fix);
gr_fix->SetMarkerStyle(21);
gr_fix->SetMarkerColor(hh->GetLineColor());
gr_fix->GetXaxis()->SetTitle("p_{t}, GeV/c");
gr_fix->SetTitle(Form("corrected phi * #sigma^{2}, %s", gtitle.Data()));
if (noSigma)
gr_fix->SetTitle(Form("corrected phi (no #sigma), %s", gtitle.Data()));
gr_fix->Draw("AP");
//---------------------
c = new TCanvas("cfinald", "data_correct", 1200, 450);
c->Divide(2, 1, 0.001, 0.001); c->Modified(); c->Draw();
c->cd(1);
gr->SetMinimum(0);
gr->SetMarkerColor(kBlack);
gr->Draw("AP");
c->cd(2);
gr_fix->SetMinimum(0);
gr_fix->SetMarkerColor(kGreen+3);
gr_fix->Draw("AP");
TString bla9 = Form("qualityonly_PID2_%s", grapht.Data());
if (noSigma) bla9 = Form("%s_noSig.C", bla9.Data());
else bla9 = Form("%s.C", bla9.Data());
// gr_fix->SaveAs(bla9.Data());
// TPad *cp = new TPad("cpf", "", 0.45,0.45,0.99,0.92);
TPad *cp = new TPad("cpf", "", 0.60,0.55,0.99,0.93);
cp->SetLogy(); cp->Draw(); cp->cd();
TGraph *cloneg = ((TGraph *)gr_fix->Clone());
cloneg->SetTitle(); cloneg->SetMarkerSize(0.8);
cloneg->Draw("AP");
// c->SaveAs(Form("%s_%s.eps", blabla.Data(), grapht.Data()));
f->Close();
}
示例9: pointsed
//.........这里部分代码省略.........
if (check == "Uncor10t2_nConstituent_") PerfTitle="uncorr pt>10 #tracks > 2 & nConstituents >1";
std::vector< TString > discriminators;
for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
discriminators.push_back( check+"disc"+taggers[itagger]+"_udsg" );
}
// discriminators.push_back( "discTC3_udsg" );
// discriminators.push_back( "discTP_udsg" );
TCanvas *cv = new TCanvas("cv","cv",900,900);
TMultiGraph *mg =new TMultiGraph();
TLegend *legend0 = new TLegend(0.68,0.12,0.88,0.32);
// inverted axis
TMultiGraph *mginv =new TMultiGraph();
TLegend *legend1 = new TLegend(0.65,0.18,0.85,0.48);
std::ofstream salida("BTagPATop3"+check+".txt");
for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
TString tag = check+"g"+taggers[itagger]+"_udsg";
TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tag);
TGraph *dgraph = (TGraph*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+discriminators[itagger]);
// TGraph *udsgvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),dgraph->GetX());
dgraph->Sort();
// udsgvsdgraph->Sort();
TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
g->Sort();
g->SetLineColor(itagger+1);
legend0 -> AddEntry(g,taggers[itagger],"l");
mg->Add(g);
//inverted axis
TGraphErrors *ginv = new TGraphErrors(agraph->GetN(),agraph->GetX(),agraph->GetY(),agraph->GetEX(),agraph->GetEY());
ginv->Sort();
ginv->SetLineColor(colorlines[itagger]);
ginv->SetMarkerStyle(8);
ginv->SetMarkerColor(colorlines[itagger]);
legend1 -> AddEntry(ginv,labeltag[itagger],"l");
mginv->Add(ginv);
std::cout << " Tagger: " << tag << std::endl;
std::cout << " Loose(10%): " << " cut > " << std::setprecision(4) << dgraph->Eval(0.1) << " b-eff = " << g->Eval(0.1) << std::endl;
std::cout << " Medium(1%): " << " cut > " << std::setprecision(4) << dgraph->Eval(0.01) << " b-eff = " << g->Eval(0.01) << std::endl;
std::cout << " Tight(0.1%) = " << " cut > " << std::setprecision(4) << dgraph->Eval(0.001) << " b-eff = " << g->Eval(0.001) << std::endl;
salida << " " << taggers[itagger] << endl;
salida << " #bin \t b-eff \t err-beff \t non-beff \t err-non-b"<< endl;
for ( size_t i=0;i< agraph->GetN();i++){
salida << " " << (i+1) << " \t "
<< agraph->GetX()[i] <<" \t "
<< agraph->GetEX()[i]<<" \t "
<< agraph->GetY()[i] <<" \t "
<< agraph->GetEY()[i]
<< " "<< endl;
}
}
// cv->SetLogx();
mg->Draw("ALP");
mg->GetYaxis()->SetTitle("b-eff");
mg->GetXaxis()->SetTitle("udsg-mistagging");
legend0 -> Draw();
cv-> SetGrid();
cv-> Print ("BTagPATop"+check+".eps");
cv-> Print ("BTagPATop"+check+".png");
cv->cd(0);
cv->SetLogx();
mg->Draw("ALP");
legend0 -> Draw();
cv-> Print ("BTagPATop1"+check+".eps");
cv-> Print ("BTagPATop1"+check+".png");
//inverted axis
TCanvas *cvinv = new TCanvas("cvinv","cvinv",700,700);
cvinv->SetLogy(0);
mginv->Draw("ALP");
mginv->GetXaxis()->SetTitle("b-eff");
mginv->GetYaxis()->SetTitle("udsg-mistagging");
legend1 -> Draw();
cvinv->SetGrid();
cvinv-> Print ("BTagPATop2"+check+".eps");
cvinv-> Print ("BTagPATop2"+check+".png");
cvinv->cd(0);
// cvinv->Update();
// cvinv->SetLogx();
// mginv->SetXmin(10^-4);
// axis range using a histogram helper
TH2D*hpx = new TH2D("hpx",PerfTitle,200,0.0,1.005,200,0.00002,1.05);
hpx->SetStats(kFALSE);
hpx->Draw();
hpx->GetXaxis()->SetTitle("b-eff");
hpx->GetYaxis()->SetTitle("udsg-mistagging");
cvinv->SetLogy(1);
mginv->Draw("AP");
legend1 -> Draw();
cvinv-> Print ("BTagPATop3"+check+".eps");
cvinv-> Print ("BTagPATop3"+check+".png");
//ed}
}
示例10: LeptonPreselectionCMG
//.........这里部分代码省略.........
if ( hardjets[j].getVarF("jn_jp") > maxJetBTag && fabs(jet.Eta()) < 2.5 )
maxJetBTag = hardjets[j].getVarF("jn_jp");
double tempDelPhiJetMet = deltaPhi(met.Phi(), jet.Phi());
if ( tempDelPhiJetMet < minDeltaPhiJetMet )
minDeltaPhiJetMet = tempDelPhiJetMet;
}
nvtx = *nvtxP;
if (isData)
ni = -1;
else
ni = *niP;
if (isSignal) {
const int nMC = ev.getSVV<int>("mcn");
const int * mcID = ev.getAVA<int>("mc_id");
int hIdx = 0;
for (; hIdx < nMC; ++hIdx)
if (fabs(mcID[hIdx]) == 25)
break;
if (hIdx == nMC)
throw string("ERROR: Higgs not found in signal sample!");
float Hpx = ev.getAVV<float>("mc_px", hIdx);
float Hpy = ev.getAVV<float>("mc_py", hIdx);
float Hpz = ev.getAVV<float>("mc_pz", hIdx);
float Hen = ev.getAVV<float>("mc_en", hIdx);
TLorentzVector higgs;
higgs.SetPxPyPzE( Hpx, Hpy, Hpz, Hen );
hmass = higgs.M();
if (higgsW) {
hweight = higgsW->Eval(hmass);
if (higgsI)
hweight *= higgsI->Eval(hmass);
} else
hweight = 1;
}
if ( opt.checkBoolOption("ADDITIONAL_LEPTON_VETO") && (type == ELE || type == MU || type == EMU) && ((nele + nmu + nsoftmu) > 2) )
continue;
if ( opt.checkBoolOption("ADDITIONAL_LEPTON_VETO") && (type == PHOT) && ((nele + nmu + nsoftmu) > 0) )
continue;
if ( opt.checkBoolOption("ZPT_CUT") && zpt < 55 )
continue;
// for different background estimation methods different window should be applied:
// * sample for photons should have 76.0 < zmass < 106.0
// * sample for non-resonant background should not have this cut applied
if ( opt.checkBoolOption("TIGHT_ZMASS_CUT") && (type == ELE || type == MU) && (zmass < 76.0 || zmass > 106.0))
continue;
if ( opt.checkBoolOption("WIDE_ZMASS_CUT") && (type == ELE || type == MU) && (zmass < 76.0 || zmass > 106.0))
continue;
if ( opt.checkBoolOption("BTAG_CUT") && ( maxJetBTag > 0.264) )
continue;
if ( opt.checkBoolOption("DPHI_CUT") && ( minDeltaPhiJetMet < 0.5) )
continue;
#ifdef PRINTEVENTS
evPrint.setJetCollection(hardjets);
evPrint.setMET(met);
evPrint.setMT(mt);
string channelType;
if (type == ELE)
channelType = "ee";
示例11: main
//.........这里部分代码省略.........
std::cerr << "... I'm minimizing ... DATA analysis" << std::endl;
std::cerr << ">>>>>>> numEvents = " << numEvents << " => " << vET_data.size() << " selected (=" << mylist->GetN() << ")" << std::endl;
numSelectedData = vET_data.size();
///===== Chi2 ====
std::cerr << " === Chi2 === " << std::endl;
minuit->SetFunction(functorChi2);
TGraph * grChi2 = new TGraph(iNoSteps);
minuit->Scan(iPar_NoBG,iNoSteps,grChi2->GetX(),grChi2->GetY(),MinScan,MaxScan);
// TGraph * grChi2 = new TGraph();
// for (int iStep = 0; iStep < iNoSteps; iStep++){
// double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
// double y = Chi2F(&x);
// grChi2->SetPoint(iStep+1,x,y);
// }
grChi2->Draw("AL");
outFile->cd();
minuit->PrintResults();
outFile->cd();
grChi2->SetTitle("grChi2");
grChi2->Write();
const double *outParametersTemp = minuit->X();
const double *errParametersTemp = minuit->Errors();
double *outParameters = new double;
double *errParameters = new double;
outParameters[0] = outParametersTemp[0];
errParameters[0] = errParametersTemp[0];
double minChi2 = grChi2->Eval(outParameters[0]);
std::cerr << " numEvents = " << numEvents << " Scale = " << outParameters[0] << " +/- " << errParameters[0] << std::endl;
///===== end Chi2 ====
///==== likelihood ====
std::cerr << " === LL === " << std::endl;
minuit->SetFunction(functorLL);
TGraph * grLL_temp = new TGraph(iNoSteps);
minuit->Scan(iPar_NoBG,iNoSteps,grLL_temp->GetX(),grLL_temp->GetY(),MinScan,MaxScan);
TGraph * grLL = new TGraph();
grLL->SetName("grLL");
int nPointLL = 0;
for (unsigned int iStep = 0; iStep < iNoSteps; iStep++){
double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
double y = LLFunc(&x);
// std::cerr << " y = " << y << std::endl;
if (y != numberDATA * numEvents) {
std::cerr << " Ok y = " << y << std::endl;
grLL->SetPoint(nPointLL,x,y);
nPointLL++;
}
}
std::cerr << " finito " << std::endl;
grLL->Draw("AL");
outFile->cd();
minuit->PrintResults();
outFile->cd();
grLL->SetTitle("grLL");
grLL->Write();
示例12: doMC_LL
///**** LL ****
void doMC_LL(){
TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
for (nIter = 0; nIter<maxIter; nIter++){
if (!(nIter%1)) std::cerr << ">>> nIter = " << nIter << " : " << maxIter << std::endl;
vET_data.clear();
outFile->cd();
TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
TH1F hDATA(nameDATA,nameDATA,numBINS,minBINS,maxBINS);
MyTreeMC->Draw(">> myListMC",(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",ScaleTrue,minET)).Data(),"entrylist");
TEntryList *myListMC = (TEntryList*)gDirectory->Get("myListMC");
MyTreeMC->SetEntryList(0);
TEntryList *listMCHere = new TEntryList("listMCHere","listMCHere");
for (int iEvt = 0; iEvt < numSelectedData; iEvt ++){
listMCHere->Enter(myListMC->GetEntry(gRandom->Uniform(0,myListMC->GetN())));
}
MyTreeMC->SetEntryList(listMCHere);
MyTreeMC->Draw(Form("(1+%f) * %s >> %s",ScaleTrue,variableName.c_str(),nameDATA.Data()));
ConvertStdVectDouble(vET_data,MyTreeMC->GetV1(),numSelectedData);
///==== likelihood ====
std::cerr << " === LL === " << std::endl;
std::cerr << " === pseudo vET_data.size() = " << vET_data.size() << std::endl;
minuit->SetFunction(functorLL);
TGraph * grLL_temp = new TGraph(iNoSteps);
minuit->Scan(iPar_NoBG,iNoSteps,grLL_temp->GetX(),grLL_temp->GetY(),MinScan,MaxScan);
TGraph * grLL = new TGraph();
int nPointLL = 0;
for (unsigned int iStep = 0; iStep < iNoSteps; iStep++){
double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
double y = LLFunc(&x);
if (y != numberDATA * numEvents) {
grLL->SetPoint(nPointLL,x,y);
nPointLL++;
}
}
grLL->Draw("AL");
outFile->cd();
minuit->PrintResults();
const double *outParametersTemp2 = minuit->X();
const double *errParametersTemp2 = minuit->Errors();
double *outParametersLL = new double;
double *errParametersLL = new double;
outParametersLL[0] = outParametersTemp2[0];
errParametersLL[0] = errParametersTemp2[0];
std::cerr << " nPointLL = " << nPointLL << std::endl;
double minLL = grLL->Eval(outParametersLL[0]);
///==== end likelihood ====
///==== Save the whole shape of LL/Chi2 ====
for (unsigned int ii=0; ii < iNoSteps; ii++){
double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
Alpha = X_ii;
Chi2 = 0;
LL = grLL->Eval(X_ii);
NewChi2 = 0;
myTreeChi2->Fill();
}
///===== Look for minima =====
double a;
double b;
double c;
double errX_low = -9999;
double errX_up = 9999;
int err_low = 0;
int err_up = 0;
for (unsigned int ii=0; ii < iNoSteps; ii++){
double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
double here = grLL->Eval(X_ii);
if (err_low == 0){
if (here < (minLL + DELTA_LL)){
errX_low = X_ii;
err_low = 1;
}
}
else if (err_up == 0 && here > (minLL + DELTA_LL) && X_ii > outParametersLL[0]){
errX_up = X_ii;
err_up = 1;
}
}
AlphaMean = outParametersLL[0];
AlphaMinus = errX_low;
AlphaPlus = errX_up;
//.........这里部分代码省略.........
示例13: bjetefficiency
void bjetefficiency() {
std::vector<string> samples_;
//samples_.push_back("WPHJET.root");
//samples_.push_back("WW.root");
//samples_.push_back("WZ.root");
//samples_.push_back("ZZ.root");
samples_.push_back("SIGNALtGC.root");
samples_.push_back("SIGNALtGu.root");
//open files
std::vector<TFile*> files;
for(unsigned int idx=0; idx<samples_.size(); ++idx){
files.push_back(new TFile(samples_[idx].c_str()));}
std::vector<string> variables2_;
/*
variables2_.push_back("btageffL");
variables2_.push_back("btageffM");
variables2_.push_back("btageffT");
variables2_.push_back("ctageffL");
variables2_.push_back("ctageffM");
variables2_.push_back("ctageffT");
*/
variables2_.push_back("btageffL");
variables2_.push_back("btageffM");
variables2_.push_back("btageffT");
std::vector<double> finaleff;
Double_t bdiscriminator[35];
for (unsigned int i=0; i<15; ++i){
bdiscriminator[i] =0.244-0.1+i*0.02; }
for (unsigned int i=0; i<10; ++i){
bdiscriminator[15+i] = 0.679-0.1+i*0.02;}
for (unsigned int i=0; i<10; ++i){
bdiscriminator[25+i] = 0.898-0.1+i*0.02;}
Double_t count[3];
count[0]=15;
count[1]=10;
count[2]=10;
for(unsigned int ii=0; ii<variables2_.size(); ++ii){
std::vector<TH1F*> hists;
for(unsigned int idx=0; idx<files.size(); ++idx){hists.push_back((TH1F*)files[idx]->Get((std::string("STEP1/").append(variables2_[ii]).c_str())));}
for(unsigned int idx=1; idx<files.size(); ++idx){
hists[idx]->Add(hists[idx-1]);}
//hists[files.size()-1]->Draw();
std::vector<double> eff;
int bin;
for ( bin=1; bin<count[ii]*2+1; ++bin){eff.push_back(hists[files.size()-1]->GetBinContent(bin));}
for (unsigned int i=0; i<count[ii]; ++i){finaleff.push_back(eff[2*i+1]/(eff[2*i]+eff[2*i+1]));
//cout<< eff[2*i+1]<<" "<<eff[2*i+1]/(eff[2*i]+eff[2*i+1])<<endl;
}
}
// TGraph *gr = new TGraph(30,bdiscriminator,finaleff);
Double_t effi[35];
for (unsigned int i=0; i<35; ++i){
effi[i] =finaleff[i]; }
TGraph *gr = new TGraph(35,effi,bdiscriminator);
TGraph *grREV = new TGraph(35,bdiscriminator,effi);
cout<<"the new working points for nominal SF"<<endl;
cout<<"loose "<<gr->Eval(1.008*grREV->Eval(0.244))<<endl;
cout<<"medium "<<gr->Eval(0.963*grREV->Eval(0.679))<<endl;
cout<<"tight "<<gr->Eval(0.947*grREV->Eval(0.898))<<endl;
cout<<"the new working points for SF UP"<<endl;
cout<<"loose "<<gr->Eval((1.008+0.023)*grREV->Eval(0.244))<<endl;
cout<<"medium "<<gr->Eval((0.963+0.020)*grREV->Eval(0.679))<<endl;
cout<<"tight "<<gr->Eval((0.947+0.025)*grREV->Eval(0.898))<<endl;
cout<<"the new working points for SF DOWN"<<endl;
cout<<"loose "<<gr->Eval((1.008-0.023)*grREV->Eval(0.244))<<endl;
cout<<"medium "<<gr->Eval((0.963-0.020)*grREV->Eval(0.679))<<endl;
cout<<"tight "<<gr->Eval((0.947-0.025)*grREV->Eval(0.898))<<endl;
// cout<<gr->Eval(1.008*grREV->Eval(0.244))<<endl;
// cout<<gr->Eval(0.966*grREV->Eval(0.679))<<endl;
// cout<<gr->Eval(0.931*grREV->Eval(0.898))<<endl;
// cout<<gr->Eval(0.0990)<<endl;
// cout<<gr->Eval(0.0142)<<endl;
// cout<<gr->Eval(0.00155)<<endl;
// cout<<gr->Eval(0.0142)<<endl;
// cout<<gr->Eval(0.0016)<<endl;
gr->SetLineColor(4);
gr->SetLineWidth(4);
gr->SetMarkerColor(4);
gr->SetMarkerStyle(21);
gr->Draw("AP");
//.........这里部分代码省略.........
示例14: main
//.........这里部分代码省略.........
if (cat==4 || cat==5 || cat==6 || cat==7 || cat==8) {
modelBuilder.addBkgPdf("Bernstein",3,Form("pol3_cat%d",cat));
modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
/*
modelBuilder.addBkgPdf("PowerLaw",1,Form("pow1_cat%d",cat));
modelBuilder.addBkgPdf("PowerLaw",3,Form("pow3_cat%d",cat));
modelBuilder.addBkgPdf("Exponential",1,Form("exp1_cat%d",cat));
modelBuilder.addBkgPdf("Exponential",3,Form("exp3_cat%d",cat));
modelBuilder.addBkgPdf("Laurent",1,Form("lau1_cat%d",cat));
modelBuilder.addBkgPdf("Laurent",3,Form("lau3_cat%d",cat));
*/
}
map<string,RooAbsPdf*> bkgPdfs = modelBuilder.getBkgPdfs();
modelBuilder.setSignalPdfFromMC(sigMC);
modelBuilder.makeSBPdfs();
map<string,RooAbsPdf*> sbPdfs = modelBuilder.getSBPdfs();
modelBuilder.fitToData(dataBinned,true,true);
modelBuilder.fitToData(dataBinned,false,true);
modelBuilder.throwToy(Form("cat%d_toy0",cat),dataBinned->sumEntries(),true,true);
// Profile this category using ProfileMultiplePdfs
ProfileMultiplePdfs profiler;
for (map<string,RooAbsPdf*>::iterator pdf=sbPdfs.begin(); pdf!=sbPdfs.end(); pdf++) {
string bkgOnlyName = pdf->first.substr(pdf->first.find("sb_")+3,string::npos);
if (bkgPdfs.find(bkgOnlyName)==bkgPdfs.end()){
cerr << "ERROR -- couldn't find bkg only pdf " << bkgOnlyName << " for SB pdf " << pdf->first << endl;
pdf->second->fitTo(*dataBinned);
exit(1);
}
int nParams = bkgPdfs[bkgOnlyName]->getVariables()->getSize()-1;
profiler.addPdf(pdf->second,2*nParams);
//profiler.addPdf(pdf->second);
cout << pdf->second->GetName() << " nParams=" << pdf->second->getVariables()->getSize() << " nBkgParams=" << nParams << endl;
}
profiler.printPdfs();
//cout << "Continue?" << endl;
//string bus; cin >> bus;
profiler.plotNominalFits(dataBinned,mass,80,Form("cat%d",cat));
pair<double,map<string,TGraph*> > minNlls = profiler.profileLikelihood(dataBinned,mass,mu,mu_low,mu_high,mu_step);
pair<double,map<string,TGraph*> > correctedNlls = profiler.computeEnvelope(minNlls,Form("cat%d",cat),2.);
minNlltrack.push_back(make_pair(correctedNlls.first,correctedNlls.second["envelope"]));
//minNlls.second.insert(pair<string,TGraph*>("envelope",envelopeNll.second));
//map<string,TGraph*> minNLLs = profiler.profileLikelihoodEnvelope(dataBinned,mu,mu_low,mu_high,mu_step);
profiler.plot(correctedNlls.second,Form("cat%d_nlls",cat));
//profiler.print(minNLLs,mu_low,mu_high,mu_step);
/*
if (minNLLs.find("envelope")==minNLLs.end()){
cerr << "ERROR -- envelope TGraph not found in minNLLs" << endl;
exit(1);
}
*/
//minNlltrack.push_back(make_pair(profiler.getGlobalMinNLL(),minNLLs["envelope"]));
}
//exit(1);
TGraph *comb = new TGraph();
for (vector<pair<double,TGraph*> >::iterator it=minNlltrack.begin(); it!=minNlltrack.end(); it++){
if (it->second->GetN()!=minNlltrack.begin()->second->GetN()){
cerr << "ERROR -- unequal number of points for TGraphs " << it->second->GetName() << " and " << minNlltrack.begin()->second->GetName() << endl;
exit(1);
}
}
for (int p=0; p<minNlltrack.begin()->second->GetN(); p++){
double x,y,sumy=0;
for (vector<pair<double,TGraph*> >::iterator it=minNlltrack.begin(); it!=minNlltrack.end(); it++){
it->second->GetPoint(p,x,y);
sumy += (y+it->first);
}
comb->SetPoint(p,x,sumy);
}
pair<double,double> globalMin = getGraphMin(comb);
for (int p=0; p<comb->GetN(); p++){
double x,y;
comb->GetPoint(p,x,y);
comb->SetPoint(p,x,y-globalMin.second);
}
vector<double> fitVal = getValsFromLikelihood(comb);
cout << "Best fit.." << endl;
cout << "\t mu = " << Form("%4.3f",fitVal[0]) << " +/- (1sig) = " << fitVal[2]-fitVal[0] << " / " << fitVal[0]-fitVal[1] << endl;
cout << "\t " << " " << " +/- (2sig) = " << fitVal[4]-fitVal[0] << " / " << fitVal[0]-fitVal[3] << endl;
cout << comb->Eval(fitVal[0]) << " " << comb->Eval(fitVal[1]) << " " << comb->Eval(fitVal[2]) << " " << comb->Eval(fitVal[3]) << " " << comb->Eval(fitVal[4]) << endl;
double quadInterpVal = ProfileMultiplePdfs::quadInterpMinimum(comb);
cout << "quadInterp: mu = " << quadInterpVal << endl;
cout << "\t " << comb->Eval(quadInterpVal) << " " << comb->Eval(quadInterpVal-0.005) << " " << comb->Eval(quadInterpVal-0.01) << " " << comb->Eval(quadInterpVal+0.005) << " " << comb->Eval(quadInterpVal+0.01) << endl;
comb->SetLineWidth(2);
TCanvas *canv = new TCanvas();
comb->Draw("ALP");
canv->Print("plots/comb.pdf");
TFile *tempOut = new TFile("tempOut.root","RECREATE");
tempOut->cd();
comb->SetName("comb");
comb->Write();
tempOut->Close();
return 0;
}
示例15: ExtractSignalYield
void ExtractSignalYield(string ModelName, string decayMode, float natural_width, bool dependsOnKtilde=0)
{
string dir = "/usr/users/dschaefer/root/results/";
string outputfile = ModelName+"_"+decayMode+"_signalYield.root";
//float natural_width =0.05;
std::stringstream s;
if(int(natural_width*100+0.5) == 10 or int(100*natural_width+0.5)==20 or int(natural_width*100+0.5)==30)
{
s << std::fixed << std::setprecision(1) << natural_width;
}
else
{
s << std::fixed << std::setprecision(2) << natural_width;
}
string swidth = s.str();
float ktilde;
string sktilde;
if(dependsOnKtilde)
{
ktilde = natural_width;
sktilde = swidth;
}
TGraph* limit = getObservedLimit(ModelName,decayMode,natural_width,dependsOnKtilde);
TGraph* limitUP = getLimitUp(ModelName,decayMode,natural_width,1,dependsOnKtilde);
TGraph* limitDOWN = getLimitDown(ModelName,decayMode,natural_width,1,dependsOnKtilde);
TGraph* limitUP2 = getLimitUp(ModelName,decayMode,natural_width,2,dependsOnKtilde);
TGraph* limitDOWN2 = getLimitDown(ModelName,decayMode,natural_width,2,dependsOnKtilde);
int number;
float* masses;
if(decayMode.find("lvjj")!=string::npos)
{
masses = new float[38];
number =38;
float array[38] = {800,900,1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000,4100,4200,4300,4400,4500};
masses = array;
}
else
{
number = 29;
masses = new float[29];
float array[29] = {1200,1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000};
masses = array;
}
TGraph* gryellow = new TGraph(2*number);
TGraph* grgreen = new TGraph(2*number);
for(int m=0;m<number;m++)
{
gryellow->SetPoint(m,masses[m],limitUP->Eval(masses[m]));
gryellow->SetPoint(number+m,masses[number-m-1],limitDOWN->Eval(masses[number-m-1]));
grgreen->SetPoint(m,masses[m],limitUP2->Eval(masses[m]));
grgreen->SetPoint(number+m,masses[number-m-1],limitDOWN2->Eval(masses[number-m-1]));
}
gryellow->SetFillColor(kYellow);
grgreen->SetFillColor(kGreen);
TCanvas* test = new TCanvas("test","test",400,400);
test->SetLogy();
limitUP->SetFillColor(kGreen+2);
limitUP->SetLineColor(kGreen+2);
limitUP->SetLineWidth(-1504);
//limitUP->Draw("ACP");
limitDOWN->SetFillColor(kBlue);
limitDOWN->SetLineColor(kGreen+2);
limitDOWN->SetLineWidth(-1504);
//limitDOWN->Draw("CPsame");
gryellow->Draw("AFsame");
grgreen->Draw("Fsame");
limit->SetLineColor(kBlue);
limit->SetLineWidth(3);
limit->Draw("Lsame");
test->Update();
float lumi = 2.10*1000; //pb^-1
if(decayMode.find("jjjj")!=string::npos){lumi = 2.6*1000;}//pb^-1
//float *xs = getTheoryCrossSection(ModelName,0.1);
//x values for which a theory prediction for production cross section exists
float m_xs[10] = {800,900,1000,1500,1800,2000,2500,3000,3500,4500};
float BR_G_to_WW[10];
float BR_G_to_ZZ[10];
float *acceptance;
float *N_expected;
//float xs = 0.01;
float *xs = new float[10];
if(dependsOnKtilde)
{
xs = getTheoryCrossSection(ModelName,ktilde);
}
//.........这里部分代码省略.........