本文整理汇总了C++中RooAbsPdf类的典型用法代码示例。如果您正苦于以下问题：C++ RooAbsPdf类的具体用法？C++ RooAbsPdf怎么用？C++ RooAbsPdf使用的例子？那么, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了RooAbsPdf类的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: mvaPUPPETEvaluation
void mvaPUPPETEvaluation() {
//output dir
TString dirname = ".";
gSystem->mkdir(dirname,true);
gSystem->cd(dirname);
//read workspace from training
TString fname;
fname = "mvaPUPPET.root";
TString infile = fname.Data();
TFile *fws = TFile::Open(infile);
RooWorkspace *ws = (RooWorkspace*)fws->Get("wereg");
//read variables from workspace
RooGBRTargetFlex *perpwidth = static_cast<RooGBRTargetFlex*>(ws->arg("perpwidth"));
RooGBRTargetFlex *perpmean = static_cast<RooGBRTargetFlex*>(ws->arg("perpmean"));
RooGBRTargetFlex *parpwidth = static_cast<RooGBRTargetFlex*>(ws->arg("parwidth"));
RooGBRTargetFlex *parmean = static_cast<RooGBRTargetFlex*>(ws->arg("parmean"));
RooRealVar *tgtvarX = ws->var("tgtX");
RooRealVar *tgtvarY = ws->var("tgtY");
RooArgList vars;
vars.add(perpwidth->FuncVars());
vars.add(perpmean->FuncVars());
vars.add(perpwidth->FuncVars());
vars.add(parmean->FuncVars());
vars.add(*tgtvarX);
vars.add(*tgtvarY);
//read testing dataset from TTree
RooRealVar weightvar("weightvar","",1.);
TChain *dtree = new TChain("tree");
dtree->Add("root://eoscms.cern.ch//store//group/dpg_ecal/alca_ecalcalib/ecalMIBI/rgerosa/PUPPETAnalysis/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8_Asympt50ns_MCRUN2_74_V9A_forMVATraining/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8/crab_20150724_111858/150724_091912/0000/output_mc_1.root/PUPPET/t");
//TFile *fdin = TFile::Open("/data/bendavid/regTreesAug1/hgg-2013Final8TeV_reg_s12-zllm50-v7n_noskim.root");
// TDirectory *ddir = (TDirectory*)fdin->FindObjectAny("");
// dtree = (TTree*)ddir->Get("t");
//selection cuts for testing
TCut selcut;
selcut = "Boson_daughter==13";
TCut selweight = "1";
TCut prescale100alt = "(evt%100==1)";
weightvar.SetTitle(selcut);
//make testing dataset
RooDataSet *hdata = RooTreeConvert::CreateDataSet("hdata",dtree,vars,weightvar);
weightvar.SetTitle(selcut);
//make reduced testing dataset for integration over conditional variables
RooDataSet *hdatasmall = RooTreeConvert::CreateDataSet("hdatasmall",dtree,vars,weightvar);
//retrieve full pdf from workspace
RooAbsPdf *sigpdfX = ws->pdf("sigpdfX");
RooAbsPdf *sigpdfY = ws->pdf("sigpdfY");
//input variable corresponding to sceta
RooRealVar *scetavar = ws->var("var_1");
//regressed output functions
RooAbsReal *perpwidthlim = ws->function("perpwidthlim");
RooAbsReal *perpmeanlim = ws->function("perpmeanlim");
RooAbsReal *parwidthlim = ws->function("parwidthlim");
RooAbsReal *parmeanlim = ws->function("parmeanlim");
//////////////////////////////////////////////////////////////////////////////////////
//
// formula for corrected recoil?
//
/*
//formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
ecorvar->setRange(0.,2.);
ecorvar->setBins(800);
//formula for raw energy/true energy (1.0/(etrue/eraw))
RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
rawvar->setRange(0.,2.);
rawvar->setBins(800);
*/
//////////////////////////////////////////////////////////////////////////////////////
//clone data and add regression outputs for plotting
RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");
RooRealVar *perpwidthvar = (RooRealVar*)hdataclone->addColumn(*perpwidthlim);
RooRealVar *perpmeanvar = (RooRealVar*)hdataclone->addColumn(*perpmeanlim);
RooRealVar *parwidthvar = (RooRealVar*)hdataclone->addColumn(*parwidthlim);
RooRealVar *parmeanvar = (RooRealVar*)hdataclone->addColumn(*parmeanlim);
//.........这里部分代码省略.........
示例2: fitbkgdataCard
void fitbkgdataCard(TString configCard="template.config",
bool dobands = true, // create baerror bands for BG models
bool dosignal = false, // plot the signal model (needs to be present)
bool blinded = true, // blind the data in the plots?
bool verbose = true ) {
gROOT->Macro("MitStyle.C");
gStyle->SetErrorX(0);
gStyle->SetOptStat(0);
gROOT->ForceStyle();
TString projectDir;
std::vector<TString> catdesc;
std::vector<TString> catnames;
std::vector<int> polorder;
double massmin = -1.;
double massmax = -1.;
double theCMenergy = -1.;
bool readStatus = readFromConfigCard( configCard,
projectDir,
catnames,
catdesc,
polorder,
massmin,
massmax,
theCMenergy
);
if( !readStatus ) {
std::cerr<<" ERROR: Could not read from card > "<<configCard.Data()<<" <."<<std::endl;
return;
}
TFile *fdata = new TFile(TString::Format("%s/CMS-HGG-data.root",projectDir.Data()),"READ");
if( !fdata ) {
std::cerr<<" ERROR: Could not open file "<<projectDir.Data()<<"/CMS-HGG-data.root."<<std::endl;
return;
}
if( !gSystem->cd(TString::Format("%s/databkg/",projectDir.Data())) ) {
std::cerr<<" ERROR: Could not change directory to "<<TString::Format("%s/databkg/",projectDir.Data()).Data()<<"."<<std::endl;
return;
}
// ----------------------------------------------------------------------
// load the input workspace....
RooWorkspace* win = (RooWorkspace*)fdata->Get("cms_hgg_workspace_data");
if( !win ) {
std::cerr<<" ERROR: Could not load workspace > cms_hgg_workspace_data < from file > "<<TString::Format("%s/CMS-HGG-data.root",projectDir.Data()).Data()<<" <."<<std::endl;
return;
}
RooRealVar *intLumi = win->var("IntLumi");
RooRealVar *hmass = win->var("CMS_hgg_mass");
if( !intLumi || !hmass ) {
std::cerr<<" ERROR: Could not load needed variables > IntLumi < or > CMS_hgg_mass < forom input workspace."<<std::endl;
return;
}
//win->Print();
hmass->setRange(massmin,massmax);
hmass->setBins(4*(int)(massmax-massmin));
hmass->SetTitle("m_{#gamma#gamma}");
hmass->setUnit("GeV");
hmass->setRange("fitrange",massmin,massmax);
hmass->setRange("blind1",100.,110.);
hmass->setRange("blind2",150.,180.);
// ----------------------------------------------------------------------
// some auxiliray vectro (don't know the meaning of all of them ... yet...
std::vector<RooAbsData*> data_vec;
std::vector<RooAbsPdf*> pdfShape_vec; // vector to store the NOT-EXTENDED PDFs (aka pdfshape)
std::vector<RooAbsPdf*> pdf_vec; // vector to store the EXTENDED PDFs
std::vector<RooAbsReal*> normu_vec; // this holds the normalization vars for each Cat (needed in bands for combined cat)
RooArgList normList; // list of range-limityed normalizations (needed for error bands on combined category)
//std::vector<RooRealVar*> coeffv;
//std::vector<RooAbsReal*> normu_vecv; // ???
// ----------------------------------------------------------------------
// define output works
RooWorkspace *wOut = new RooWorkspace("wbkg","wbkg") ;
// util;ities for the combined fit
RooCategory finalcat ("finalcat", "finalcat") ;
RooSimultaneous fullbkgpdf("fullbkgpdf","fullbkgpdf",finalcat);
RooDataSet datacomb ("datacomb", "datacomb", RooArgList(*hmass,finalcat)) ;
RooDataSet *datacombcat = new RooDataSet("data_combcat","",RooArgList(*hmass)) ;
// add the 'combcat' to the list...if more than one cat
if( catnames.size() > 1 ) {
//.........这里部分代码省略.........
示例3: eregtesting_13TeV_Pi0_lessP2
//.........这里部分代码省略.........
TCut prescale10alt = "(Entry$%10==1)";
TCut prescale25 = "(Entry$%25==0)";
TCut prescale100 = "(Entry$%100==0)";
TCut prescale1000 = "(Entry$%1000==0)";
TCut evenevents = "(Entry$%2==0)";
TCut oddevents = "(Entry$%2==1)";
TCut prescale100alt = "(Entry$%100==1)";
TCut prescale1000alt = "(Entry$%1000==1)";
TCut prescale50alt = "(Entry$%50==1)";
TCut Events3_4 = "(Entry$%4==3)";
TCut Events1_4 = "(Entry$%4==1)";
TCut Events2_4 = "(Entry$%4==2)";
TCut Events0_4 = "(Entry$%4==0)";
TCut Events01_4 = "(Entry$%4<2)";
TCut Events23_4 = "(Entry$%4>1)";
if (doele)
weightvar.SetTitle(prescale100alt*selcut);
else
weightvar.SetTitle(Events01_4*selcut);
//make testing dataset
RooDataSet *hdata = RooTreeConvert::CreateDataSet("hdata",dtree,vars,weightvar);
if (doele)
weightvar.SetTitle(prescale1000alt*selcut);
else
weightvar.SetTitle(prescale10alt*selcut);
//make reduced testing dataset for integration over conditional variables
RooDataSet *hdatasmall = RooTreeConvert::CreateDataSet("hdatasmall",dtree,vars,weightvar);
//retrieve full pdf from workspace
RooAbsPdf *sigpdf = ws->pdf("sigpdf");
//input variable corresponding to sceta
//// RooRealVar *scetavar = ws->var("var_1");
//// RooRealVar *scphivar = ws->var("var_2");
//regressed output functions
RooAbsReal *sigmeanlim = ws->function("sigmeanlim");
RooAbsReal *sigwidthlim = ws->function("sigwidthlim");
RooAbsReal *signlim = ws->function("signlim");
RooAbsReal *sign2lim = ws->function("sign2lim");
//formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
ecorvar->setRange(0.,2.);
ecorvar->setBins(800);
//formula for raw energy/true energy (1.0/(etrue/eraw))
RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
rawvar->setRange(0.,2.);
rawvar->setBins(800);
//clone data and add regression outputs for plotting
RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");
RooRealVar *meanvar = (RooRealVar*)hdataclone->addColumn(*sigmeanlim);
RooRealVar *widthvar = (RooRealVar*)hdataclone->addColumn(*sigwidthlim);
RooRealVar *nvar = (RooRealVar*)hdataclone->addColumn(*signlim);
RooRealVar *n2var = (RooRealVar*)hdataclone->addColumn(*sign2lim);
示例4: Tbroomfit
//.........这里部分代码省略.........
RooRealVar tail1("tail1", "tail", tailstart, tailmin, tailmax );msat[3][0]=&tail1;
RooRealVar tail2("tail2", "tail", tailstart, tailmin, tailmax );msat[3][1]=&tail2;
RooRealVar tail3("tail3", "tail", tailstart, tailmin, tailmax );msat[3][2]=&tail3;
RooRealVar tail4("tail4", "tail", tailstart, tailmin, tailmax );msat[3][3]=&tail4;
RooRealVar tail5("tail5", "tail", tailstart, tailmin, tailmax );msat[3][4]=&tail5;
RooRealVar tail6("tail6", "tail", tailstart, tailmin, tailmax );msat[3][5]=&tail6;
// for CBShape
RooRealVar nalpha1("nalpha1", "nalpha", 1.3, 0, 100 );msat[4][0]=&nalpha1;
RooRealVar n1("n1", "n", 5.1, 0, 100 ); msat[5][0]=&n1;
/*
* initial values for peak positions................
*/
if (npeaks>=1) {mean1=peak[0];sigma1=sigm[0];}
if (npeaks>=2) {mean2=peak[1];sigma2=sigm[1];}
if (npeaks>=3) {mean3=peak[2];sigma3=sigm[2];}
if (npeaks>=4) {mean4=peak[3];sigma4=sigm[3];}
if (npeaks>=5) {mean5=peak[4];sigma5=sigm[4];}
if (npeaks>=6) {mean6=peak[5];sigma6=sigm[5];}
/*
* RooAbsPdf -> RooGaussian
* RooNovosibirsk
* RooLandau
*/
RooAbsPdf *pk[6]; // MAXIMUM PEAKS ==5 6 NOW!!
RooAbsPdf *pk_dicto[14][6]; // ALL DICTIONARY OF PEAKS..........
// Abstract Class.... carrefuly
RooGaussian gauss1("gauss1","gauss(x,mean,sigma)", x, mean1, sigma1);pk_dicto[0][0]=&gauss1;
RooGaussian gauss2("gauss2","gauss(x,mean,sigma)", x, mean2, sigma2);pk_dicto[0][1]=&gauss2;
RooGaussian gauss3("gauss3","gauss(x,mean,sigma)", x, mean3, sigma3);pk_dicto[0][2]=&gauss3;
RooGaussian gauss4("gauss4","gauss(x,mean,sigma)", x, mean4, sigma4);pk_dicto[0][3]=&gauss4;
RooGaussian gauss5("gauss5","gauss(x,mean,sigma)", x, mean5, sigma5);pk_dicto[0][4]=&gauss5;
RooGaussian gauss6("gauss6","gauss(x,mean,sigma)", x, mean6, sigma6);pk_dicto[0][5]=&gauss6;
RooNovosibirsk ns1("ns1","novosib(x,mean,sigma,tail)", x, mean1,sigma1, tail1 );pk_dicto[1][0]=&ns1;
RooNovosibirsk ns2("ns2","novosib(x,mean,sigma,tail)", x, mean2,sigma2, tail2 );pk_dicto[1][1]=&ns2;
RooNovosibirsk ns3("ns3","novosib(x,mean,sigma,tail)", x, mean3,sigma3, tail3 );pk_dicto[1][2]=&ns3;
RooNovosibirsk ns4("ns4","novosib(x,mean,sigma,tail)", x, mean4,sigma4, tail4 );pk_dicto[1][3]=&ns4;
RooNovosibirsk ns5("ns5","novosib(x,mean,sigma,tail)", x, mean5,sigma5, tail5 );pk_dicto[1][4]=&ns5;
// BreitWiegner is Lorentzian...?
RooBreitWigner bw1("bw1","BreitWigner(x,mean,sigma)", x, mean1, sigma1 );pk_dicto[2][0]=&bw1;
RooBreitWigner bw2("bw2","BreitWigner(x,mean,sigma)", x, mean2, sigma2 );pk_dicto[2][1]=&bw2;
RooBreitWigner bw3("bw3","BreitWigner(x,mean,sigma)", x, mean3, sigma3 );pk_dicto[2][2]=&bw3;
RooBreitWigner bw4("bw4","BreitWigner(x,mean,sigma)", x, mean4, sigma4 );pk_dicto[2][3]=&bw4;
RooBreitWigner bw5("bw5","BreitWigner(x,mean,sigma)", x, mean5, sigma5 );pk_dicto[2][4]=&bw5;
RooCBShape cb1("cb1","CBShape(x,mean,sigma)", x, mean1, sigma1, nalpha1, n1 );pk_dicto[3][0]=&cb1;
RooCBShape cb2("cb2","CBShape(x,mean,sigma)", x, mean2, sigma2, nalpha1, n1 );pk_dicto[3][1]=&cb2;
RooCBShape cb3("cb3","CBShape(x,mean,sigma)", x, mean3, sigma3, nalpha1, n1 );pk_dicto[3][2]=&cb3;
RooCBShape cb4("cb4","CBShape(x,mean,sigma)", x, mean4, sigma4, nalpha1, n1 );pk_dicto[3][3]=&cb4;
RooCBShape cb5("cb5","CBShape(x,mean,sigma)", x, mean5, sigma5, nalpha1, n1 );pk_dicto[3][4]=&cb5;
RooCBShape cb6("cb6","CBShape(x,mean,sigma)", x, mean6, sigma6, nalpha1, n1 );pk_dicto[3][5]=&cb6;
示例5: fit2015
void fit2015(
TString FileName ="/afs/cern.ch/user/a/anstahll/work/public/ExpressStream2015/ppData/OniaTree_262163_262328.root",
int oniamode = 2, // oniamode-> 3: Z, 2: Upsilon and 1: J/Psi
bool isData = true, // isData = false for MC, true for Data
bool isPbPb = false, // isPbPb = false for pp, true for PbPb
bool doFit = false ,
bool inExcStat = true // if inExcStat is true, then the excited states are fitted
) {
InputOpt opt;
SetOptions(&opt, isData, isPbPb, oniamode,inExcStat);
if (isPbPb) {
FileName = "/afs/cern.ch/user/a/anstahll/work/public/ExpressStream2015/PbPbData/OniaTree_262548_262893.root";
} else {
FileName = "/afs/cern.ch/user/a/anstahll/work/public/ExpressStream2015/ppData/OniaTree_262163_262328.root";
}
int nbins = 1; //ceil((opt.dMuon->M->Max - opt.dMuon->M->Min)/binw);
if (oniamode==1){
nbins = 140;
} else if (oniamode==2) {
nbins = 70;
} else if (oniamode==3) {
nbins = 40;
}
RooWorkspace myws;
TH1F* hDataOS = new TH1F("hDataOS","hDataOS", nbins, opt.dMuon.M.Min, opt.dMuon.M.Max);
makeWorkspace2015(myws, FileName, opt, hDataOS);
RooRealVar* mass = (RooRealVar*) myws.var("invariantMass");
RooDataSet* dataOS_fit = (RooDataSet*) myws.data("dataOS");
RooDataSet* dataSS_fit = (RooDataSet*) myws.data("dataSS");
RooAbsPdf* pdf = NULL;
if (oniamode==3) { doFit=false; }
if (doFit) {
int sigModel=0, bkgModel=0;
if (isData) {
if (oniamode==1){
sigModel = inExcStat ? 2 : 3;
bkgModel = 1;
} else {
sigModel = inExcStat ? 1 : 3; // gaussian
bkgModel = 2;
}
} else {
if (oniamode==1){
sigModel = inExcStat ? 2 : 3; // gaussian
bkgModel = 2;
} else {
sigModel = inExcStat ? 2 : 3; // gaussian
bkgModel = 3;
}
}
if (opt.oniaMode==1) buildModelJpsi2015(myws, sigModel, bkgModel,inExcStat);
else if (opt.oniaMode==2) buildModelUpsi2015(myws, sigModel, bkgModel,inExcStat);
pdf =(RooAbsPdf*) myws.pdf("pdf");
RooFitResult* fitObject = pdf->fitTo(*dataOS_fit,Save(),Hesse(kTRUE),Extended(kTRUE)); // Fit
}
RooPlot* frame = mass->frame(Bins(nbins),Range(opt.dMuon.M.Min, opt.dMuon.M.Max));
RooPlot* frame2 = NULL;
dataSS_fit->plotOn(frame, Name("dataSS_FIT"), MarkerColor(kRed), LineColor(kRed), MarkerSize(1.2));
dataOS_fit->plotOn(frame, Name("dataOS_FIT"), MarkerColor(kBlue), LineColor(kBlue), MarkerSize(1.2));
if (doFit) {
pdf->plotOn(frame,Name("thePdf"),Normalization(dataOS_fit->sumEntries(),RooAbsReal::NumEvent));
RooHist *hpull = frame -> pullHist(0,0,true);
hpull -> SetName("hpull");
frame2 = mass->frame(Title("Pull Distribution"),Bins(nbins),Range(opt.dMuon.M.Min,opt.dMuon.M.Max));
frame2 -> addPlotable(hpull,"PX");
}
drawPlot(frame,frame2, pdf, opt, doFit,inExcStat);
TString OutputFileName = "";
if (isPbPb) {
FileName = "/afs/cern.ch/user/a/anstahll/work/public/ExpressStream2015/PbPbData/OniaTree_262548_262893.root";
opt.RunNb.Start=262548;
opt.RunNb.End=262893;
if (oniamode==1) {OutputFileName = (TString)("JPSIPbPbDataset.root");}
if (oniamode==2) {OutputFileName = (TString)("YPbPbDataset.root");}
if (oniamode==3) {OutputFileName = (TString)("ZPbPbDataset.root");}
} else {
FileName = "/afs/cern.ch/user/a/anstahll/work/public/ExpressStream2015/ppData/OniaTree_262163_262328.root";
opt.RunNb.Start=262163;
opt.RunNb.End=262328;
if (oniamode==1) {OutputFileName = (TString)("JPSIppDataset.root");}
if (oniamode==2) {OutputFileName = (TString)("YppDataset.root");}
if (oniamode==3) {OutputFileName = (TString)("ZppDataset.root");}
}
TFile* oFile = new TFile(OutputFileName,"RECREATE");
oFile->cd();
hDataOS->Write("hDataOS");
dataOS_fit->Write("dataOS_FIT");
//.........这里部分代码省略.........
示例6: LeptonPreselectionCMG
//.........这里部分代码省略.........
smallTree->Branch( "L2ETA", &l2eta, "L2ETA/D" );
smallTree->Branch( "L2PHI", &l2phi, "L2PHI/D" );
smallTree->Branch( "ZMASS", &zmass, "ZMASS/D" );
smallTree->Branch( "ZPT", &zpt, "ZPT/D" );
smallTree->Branch( "ZETA", &zeta, "ZETA/D" );
smallTree->Branch( "MT", &mt, "MT/D" );
smallTree->Branch( "NSOFTJET", &nsoftjet, "NSOFTJET/I" );
smallTree->Branch( "NHARDJET", &nhardjet, "NHARDJET/I" );
smallTree->Branch( "MAXJETBTAG", &maxJetBTag, "MAXJETBTAG/D" );
smallTree->Branch( "MINDPJETMET", &minDeltaPhiJetMet, "MINDPJETMET/D" );
smallTree->Branch( "DETAJJ", &detajj, "DETAJJ/D" );
smallTree->Branch( "MJJ", &mjj, "MJJ/D" );
smallTree->Branch( "NVTX", &nvtx, "NVTX/I" );
smallTree->Branch( "nInter" , &ni, "nInter/I" );
smallTree->Branch( "CATEGORY", &category, "CATEGORY/I" );
smallTree->Branch( "Weight" , &weight, "Weight/D" );
smallTree->Branch( "HMASS", &hmass, "HMASS/D" );
smallTree->Branch( "HWEIGHT", &hweight, "HWEIGHT/D" );
bool isData = opt.checkBoolOption("DATA");
unsigned long nentries = tree->GetEntries();
RooDataSet * events = nullptr;
PhotonPrescale photonPrescales;
vector<int> thresholds;
if (type == PHOT) {
if (w == nullptr)
throw string("ERROR: No mass peak pdf!");
RooRealVar * zmass = w->var("mass");
zmass->setRange(76.0, 106.0);
RooAbsPdf * pdf = w->pdf("massPDF");
events = pdf->generate(*zmass, nentries);
photonPrescales.addTrigger("HLT_Photon36_R9Id90_HE10_Iso40_EBOnly", 36, 3, 7);
photonPrescales.addTrigger("HLT_Photon50_R9Id90_HE10_Iso40_EBOnly", 50, 5, 8);
photonPrescales.addTrigger("HLT_Photon75_R9Id90_HE10_Iso40_EBOnly", 75, 7, 9);
photonPrescales.addTrigger("HLT_Photon90_R9Id90_HE10_Iso40_EBOnly", 90, 10, 10);
}
TH1D ptSpectrum("ptSpectrum", "ptSpectrum", 200, 55, 755);
ptSpectrum.Sumw2();
unordered_set<EventAdr> eventsSet;
for ( unsigned long iEvent = 0; iEvent < nentries; iEvent++ ) {
// if (iEvent < 6060000)
// continue;
if ( iEvent % 10000 == 0) {
cout << string(40, '\b');
cout << setw(10) << iEvent << " / " << setw(10) << nentries << " done ..." << std::flush;
}
tree->GetEntry( iEvent );
run = -999;
lumi = -999;
event = -999;
pfmet = -999;
nele = -999;
nmu = -999;
nsoftmu = -999;
l1pt = -999;
l1eta = -999;示例7: Error
//.........这里部分代码省略.........
bModel = (ModelConfig*) sbModel->Clone();
bModel->SetName(TString(modelSBName)+TString("_with_poi_0"));
RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
if (!var) return 0;
double oldval = var->getVal();
var->setVal(0);
bModel->SetSnapshot( RooArgSet(*var) );
var->setVal(oldval);
}
else {
if (!bModel->GetSnapshot() ) {
Info("StandardHypoTestInvDemo","Model %s has no snapshot - make one using model poi and 0 values ",modelBName);
RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
if (var) {
double oldval = var->getVal();
var->setVal(0);
bModel->SetSnapshot( RooArgSet(*var) );
var->setVal(oldval);
}
else {
Error("StandardHypoTestInvDemo","Model %s has no valid poi",modelBName);
return 0;
}
}
}
// check model has global observables when there are nuisance pdf
// for the hybrid case the globobs are not needed
if (type != 1 ) {
bool hasNuisParam = (sbModel->GetNuisanceParameters() && sbModel->GetNuisanceParameters()->getSize() > 0);
bool hasGlobalObs = (sbModel->GetGlobalObservables() && sbModel->GetGlobalObservables()->getSize() > 0);
if (hasNuisParam && !hasGlobalObs ) {
// try to see if model has nuisance parameters first
RooAbsPdf * constrPdf = RooStats::MakeNuisancePdf(*sbModel,"nuisanceConstraintPdf_sbmodel");
if (constrPdf) {
Warning("StandardHypoTestInvDemo","Model %s has nuisance parameters but no global observables associated",sbModel->GetName());
Warning("StandardHypoTestInvDemo","\tThe effect of the nuisance parameters will not be treated correctly ");
}
}
}
// save all initial parameters of the model including the global observables
RooArgSet initialParameters;
RooArgSet * allParams = sbModel->GetPdf()->getParameters(*data);
allParams->snapshot(initialParameters);
delete allParams;
// run first a data fit
const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
RooRealVar *poi = (RooRealVar*)poiSet->first();
std::cout << "StandardHypoTestInvDemo : POI initial value: " << poi->GetName() << " = " << poi->getVal() << std::endl;
// fit the data first (need to use constraint )
TStopwatch tw;
bool doFit = initialFit;
if (testStatType == 0 && initialFit == -1) doFit = false; // case of LEP test statistic
if (type == 3 && initialFit == -1) doFit = false; // case of Asymptoticcalculator with nominal Asimov
double poihat = 0;
if (minimizerType.size()==0) minimizerType = ROOT::Math::MinimizerOptions::DefaultMinimizerType();
else
ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizerType.c_str());
示例8: MakePlots
void MakePlots(RooWorkspace* ws){
std::cout << "make plots" << std::endl;
RooAbsPdf* model = ws->pdf("model");
RooAbsPdf* model_dY = ws->pdf("model_dY");
RooRealVar* nsig = ws->var("nsig");
RooRealVar* nsigDPS = ws->var("nsigDPS");
RooRealVar* nsigSPS = ws->var("nsigSPS");
RooRealVar* nBbkg = ws->var("nBbkg");
RooRealVar* nbkg = ws->var("nbkg");
RooRealVar* nbkg2 = ws->var("nbkg2");
RooRealVar* FourMu_Mass = ws->var("FourMu_Mass");
RooRealVar* Psi1_Mass = ws->var("Psi1_Mass");
RooRealVar* Psi2_Mass = ws->var("Psi2_Mass");
RooRealVar* Psi1_CTxy = ws->var("Psi1_CTxy");
RooRealVar* Psi2_CTxy = ws->var("Psi2_CTxy");
RooRealVar* Psi1To2_dY = ws->var("Psi1To2_dY");
RooRealVar* Psi1To2Significance = ws->var("Psi1To2Significance");
// note, we get the dataset with sWeights
RooDataSet* data = (RooDataSet*) ws->data("dataWithSWeights");
model->fitTo(*data, Extended() );
// make TTree with efficiency variation info
//TFile *fFile = new TFile("Fit_Results_Eff_Cut.root","recreate");
setTDRStyle();
// make our canvas
TCanvas* cmass1 = new TCanvas("sPlotMass1","sPlotMass1", 600, 600);
cmass1->cd();
cmass1->SetFillColor(kWhite);
RooPlot* Mass1Plot = Psi1_Mass->frame(20);
data->plotOn(Mass1Plot,Name("data"), DataError(RooAbsData::SumW2));
model->plotOn(Mass1Plot,Name("all"));
model->plotOn(Mass1Plot,Name("sig"),RooFit::Components("Sig,sig_*"),RooFit::LineColor(kRed), RooFit::LineStyle(2));
model->plotOn(Mass1Plot,Name("bkg"),RooFit::Components("bkg_model"),RooFit::LineColor(kGreen), RooFit::LineStyle(3));
model->plotOn(Mass1Plot,Name("bkg2"),RooFit::Components("bkg2_model"),RooFit::LineColor(kBlack), RooFit::LineStyle(4));
model->plotOn(Mass1Plot,Name("Bbkg"),RooFit::Components("BBkg,Bbkg_*"),RooFit::LineColor(6), RooFit::LineStyle(7));
Mass1Plot->SetTitle("");
Mass1Plot->SetXTitle("#mu^{+}#mu^{-} 1 Invariant Mass (GeV/c^{2})");
Mass1Plot->SetYTitle("Events / 0.025 GeV/c^{2}");
Mass1Plot->SetLabelOffset(0.012);
//Mass1Plot->SetTitleOffset(0.95);
Mass1Plot->Draw();
//cmass1->SaveAs("pic/Psi1_mass.pdf");
//cmass1->Close();
TCanvas* cmass2 = new TCanvas("sPlotMass2","sPlotMass2", 600, 600);
cmass2->cd();
cmass2->SetFillColor(kWhite);
RooPlot* Mass2Plot = Psi2_Mass->frame(20);
data->plotOn(Mass2Plot,Name("data"), DataError(RooAbsData::SumW2));
model->plotOn(Mass2Plot,Name("all"));
model->plotOn(Mass2Plot,Name("sig"),RooFit::Components("Sig,sig_*"),RooFit::LineColor(kRed), RooFit::LineStyle(2));
model->plotOn(Mass2Plot,Name("bkg"),RooFit::Components("bkg_model"),RooFit::LineColor(kGreen), RooFit::LineStyle(3));
model->plotOn(Mass2Plot,Name("bkg2"),RooFit::Components("bkg2_model"),RooFit::LineColor(kBlack), RooFit::LineStyle(4));
model->plotOn(Mass2Plot,Name("Bbkg"),RooFit::Components("BBkg,Bbkg_*"),RooFit::LineColor(6), RooFit::LineStyle(7));
Mass2Plot->SetTitle("");
Mass2Plot->SetXTitle("#mu^{+}#mu^{-} 2 Invariant Mass (GeV/c^{2})");
Mass2Plot->SetYTitle("Events / 0.025 GeV/c^{2}");
Mass2Plot->SetLabelOffset(0.012);
//Mass2Plot->SetTitleOffset(0.95);
Mass2Plot->Draw();
//cmass2->SaveAs("pic/Psi2_mass.pdf");
//cmass2->Close();
TCanvas* cctxy1 = new TCanvas("sPlotCTxy1","sPlotCTxy1", 600, 600);
cctxy1->cd();
cctxy1->SetFillColor(kWhite);
RooPlot* CTxy1Plot = Psi1_CTxy->frame(30);
data->plotOn(CTxy1Plot,Name("data"), DataError(RooAbsData::SumW2));
model->plotOn(CTxy1Plot,Name("all"));
model->plotOn(CTxy1Plot,Name("sig"),RooFit::Components("Sig,sig_*"),RooFit::LineColor(kRed), RooFit::LineStyle(2));
model->plotOn(CTxy1Plot,Name("bkg"),RooFit::Components("bkg_model"),RooFit::LineColor(kGreen), RooFit::LineStyle(3));
model->plotOn(CTxy1Plot,Name("bkg2"),RooFit::Components("bkg2_model"),RooFit::LineColor(kBlack), RooFit::LineStyle(4));
model->plotOn(CTxy1Plot,Name("Bbkg"),RooFit::Components("BBkg,Bbkg_*"),RooFit::LineColor(6), RooFit::LineStyle(7));
CTxy1Plot->SetTitle("");
CTxy1Plot->SetXTitle("J/#psi^{1} ct_{xy} (cm)");
CTxy1Plot->SetYTitle("Events / 0.005 cm");
CTxy1Plot->SetMaximum(2000);
CTxy1Plot->SetMinimum(0.1);
CTxy1Plot->Draw();
cctxy1->SetLogy();
//cctxy1->SaveAs("pic/Psi1_CTxy.pdf");
//cctxy1->Close();
TCanvas* csig = new TCanvas("sPlotSig","sPlotSig", 600, 600);
csig->cd();
RooPlot* SigPlot = Psi1To2Significance->frame(20);
data->plotOn(SigPlot,Name("data"), DataError(RooAbsData::SumW2));
model->plotOn(SigPlot,Name("all"));
model->plotOn(SigPlot,Name("sig"),RooFit::Components("Sig,sig_*"),RooFit::LineColor(kRed), RooFit::LineStyle(2));
model->plotOn(SigPlot,Name("bkg"),RooFit::Components("bkg_model"),RooFit::LineColor(kGreen), RooFit::LineStyle(3));
model->plotOn(SigPlot,Name("bkg2"),RooFit::Components("bkg2_model"),RooFit::LineColor(kBlack), RooFit::LineStyle(4));
model->plotOn(SigPlot,Name("Bbkg"),RooFit::Components("BBkg,Bbkg_*"),RooFit::LineColor(6), RooFit::LineStyle(7));
SigPlot->SetTitle("");
SigPlot->SetYTitle("Events / 0.4");
SigPlot->SetXTitle("J/#psi Distance Significance");
SigPlot->Draw();
//.........这里部分代码省略.........
示例9: ComputeTestStat
float ComputeTestStat(TString wsfile, double mu_susy_sig_val) {
gROOT->Reset();
TFile* wstf = new TFile( wsfile ) ;
RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );
ws->Print() ;
ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
modelConfig->Print() ;
RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;
rds->Print() ;
rds->printMultiline(cout, 1, kTRUE, "") ;
RooAbsPdf* likelihood = modelConfig->GetPdf() ;
RooRealVar* rrv_mu_susy_sig = ws->var("mu_susy_all0lep") ;
if ( rrv_mu_susy_sig == 0x0 ) {
printf("\n\n\n *** can't find mu_susy_all0lep in workspace. Quitting.\n\n\n") ;
return ;
} else {
printf(" current value is : %8.3f\n", rrv_mu_susy_sig->getVal() ) ; cout << flush ;
rrv_mu_susy_sig->setConstant(kFALSE) ;
}
/*
// check the impact of varying the qcd normalization:
RooRealVar *rrv_qcd_0lepLDP_ratioH1 = ws->var("qcd_0lepLDP_ratio_H1");
RooRealVar *rrv_qcd_0lepLDP_ratioH2 = ws->var("qcd_0lepLDP_ratio_H2");
RooRealVar *rrv_qcd_0lepLDP_ratioH3 = ws->var("qcd_0lepLDP_ratio_H3");
rrv_qcd_0lepLDP_ratioH1->setVal(0.3);
rrv_qcd_0lepLDP_ratioH2->setVal(0.3);
rrv_qcd_0lepLDP_ratioH3->setVal(0.3);
rrv_qcd_0lepLDP_ratioH1->setConstant(kTRUE);
rrv_qcd_0lepLDP_ratioH2->setConstant(kTRUE);
rrv_qcd_0lepLDP_ratioH3->setConstant(kTRUE);
*/
printf("\n\n\n ===== Doing a fit with SUSY component floating ====================\n\n") ;
RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
double logLikelihoodSusyFloat = fitResult->minNll() ;
double logLikelihoodSusyFixed(0.) ;
double testStatVal(-1.) ;
if ( mu_susy_sig_val >= 0. ) {
printf("\n\n\n ===== Doing a fit with SUSY fixed ====================\n\n") ;
printf(" fixing mu_susy_sig to %8.2f.\n", mu_susy_sig_val ) ;
rrv_mu_susy_sig->setVal( mu_susy_sig_val ) ;
rrv_mu_susy_sig->setConstant(kTRUE) ;
fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
logLikelihoodSusyFixed = fitResult->minNll() ;
testStatVal = 2.*(logLikelihoodSusyFixed - logLikelihoodSusyFloat) ;
printf("\n\n\n ======= test statistic : -2 * ln (L_fixed / ln L_max) = %8.3f\n\n\n", testStatVal ) ;
}
return testStatVal ;
}示例10: drawSidebandsWithCurve
void drawSidebandsWithCurve( DrawBase* db, const std::string& data_dataset, const std::string& PUType, const std::string& data_mc, const std::string& flags, int nbtags, const std::string& leptType, std::string suffix ) {
if( suffix!="" ) suffix = "_" + suffix;
TH1F::AddDirectory(kTRUE);
// get histograms:
std::vector< TH1D* > lastHistos_data = db->get_lastHistos_data();
TH1D* h1_data = new TH1D(*(lastHistos_data[0]));
float xMin = (db->get_xAxisMin()!=9999.) ? db->get_xAxisMin() : h1_data->GetXaxis()->GetXmin();
float xMax = (db->get_xAxisMax()!=9999.) ? db->get_xAxisMax() : h1_data->GetXaxis()->GetXmax();
int nBins = (int)(xMax-xMin)/h1_data->GetXaxis()->GetBinWidth(1);
std::string sherpa_suffix = (use_sherpa) ? "_sherpa" : "";
// open fit results file:
char fitResultsFileName[200];
sprintf( fitResultsFileName, "fitResultsFile_%s_%dbtag_ALL_PU%s_fit%s%s%s.root", data_dataset.c_str(), nbtags, PUType.c_str(), data_mc.c_str(), sherpa_suffix.c_str(), flags.c_str() );
TFile* fitResultsFile = TFile::Open(fitResultsFileName);
// get bg workspace:
char workspaceName[200];
sprintf( workspaceName, "fitWorkspace_%dbtag", nbtags );
RooWorkspace* bgws = (RooWorkspace*)fitResultsFile->Get(workspaceName);
char fitResultsName[200];
sprintf( fitResultsName, "fitResults_%dbtag_decorr", nbtags );
RooFitResult* bgfr = (RooFitResult*)fitResultsFile->Get(fitResultsName);
// get mZZ variable:
RooRealVar* CMS_hzz2l2q_mZZ = (RooRealVar*)bgws->var("CMS_hzz2l2q_mZZ");
// get bg shape:
RooAbsPdf* background = (RooAbsPdf*)bgws->pdf("background_decorr");
// get data sidebands:
TTree* tree_sidebands = (TTree*)fitResultsFile->Get("sidebandsDATA_alpha");
TH1D* h1_dataSidebands = new TH1D("dataSidebands", "", nBins, xMin, xMax);
h1_dataSidebands->Sumw2();
char selection[900];
if( leptType=="ALL" )
sprintf( selection, "eventWeight_alpha*(((mZjj>60.&&mZjj<75.)||(mZjj>105.&&mZjj<130.)) && nBTags==%d && CMS_hzz2l2q_mZZ>183. && CMS_hzz2l2q_mZZ<800.)", nbtags);
else {
int leptType_int = SidebandFitter::convert_leptType(leptType);
sprintf( selection, "eventWeight_alpha*(((mZjj>60.&&mZjj<75.)||(mZjj>105.&&mZjj<130.)) && nBTags==%d && leptType==%d && CMS_hzz2l2q_mZZ>183. && CMS_hzz2l2q_mZZ<800.)", nbtags, leptType_int);
}
tree_sidebands->Project("dataSidebands", "CMS_hzz2l2q_mZZ", selection);
// create data graph (poisson asymm errors):
TGraphAsymmErrors* graph_data_poisson = new TGraphAsymmErrors(0);
graph_data_poisson = fitTools::getGraphPoissonErrors(h1_dataSidebands);
graph_data_poisson->SetMarkerStyle(20);
std::string leptType_cut;
if( leptType=="MU" ) leptType_cut = "&& leptType==0";
if( leptType=="ELE" ) leptType_cut = "&& leptType==1";
SidebandFitter* sf = new SidebandFitter( data_dataset, PUType, data_mc, flags );
// get bg normalization:
//float expBkg = sf->get_backgroundNormalization( nbtags, leptType, "DATA" );
std::pair<Double_t, Double_t> bgNormAndError = sf->get_backgroundNormalizationAndError( nbtags, leptType, "DATA" );
float expBkg = bgNormAndError.first;
float expBkg_err = bgNormAndError.second;
//TTree* treeSidebandsDATA_alphaCorr = (TTree*)fitResultsFile->Get("sidebandsDATA_alpha");
//TH1D* h1_mZZ_sidebands_alpha = new TH1D("mZZ_sidebands_alpha", "", 65, xMin, xMax);
//char sidebandsCut_alpha[500];
//sprintf(sidebandsCut_alpha, "eventWeight_alpha*(isSidebands && nBTags==%d %s)", nbtags, leptType_cut.c_str());
//treeSidebandsDATA_alphaCorr->Project("mZZ_sidebands_alpha", "CMS_hzz2l2q_mZZ", sidebandsCut_alpha);
//float expBkg = h1_mZZ_sidebands_alpha->Integral();
RooPlot *plot_MCbkg = CMS_hzz2l2q_mZZ->frame(xMin,xMax,nBins);
background->plotOn(plot_MCbkg,RooFit::Normalization(expBkg));
if( !QUICK_ ) {
background->plotOn(plot_MCbkg,RooFit::VisualizeError(*bgfr,2.0,kFALSE),RooFit::FillColor(kYellow), RooFit::Normalization(expBkg));
background->plotOn(plot_MCbkg,RooFit::VisualizeError(*bgfr,1.0,kFALSE),RooFit::FillColor(kGreen), RooFit::Normalization(expBkg));
background->plotOn(plot_MCbkg,RooFit::Normalization(expBkg));
}
TF1* f1_bgForLegend = new TF1("bgForLegend", "[0]");
f1_bgForLegend->SetLineColor(kBlue);
f1_bgForLegend->SetLineWidth(3);
TH2D* h2_axes = new TH2D("axes", "", 10, xMin, xMax, 10, 0., 1.3*h1_data->GetMaximum());
//.........这里部分代码省略.........
示例11: DoSPlot
//____________________________________
void DoSPlot(RooWorkspace* ws){
std::cout << "Calculate sWeights" << std::endl;
RooAbsPdf* model = ws->pdf("model");
RooRealVar* nsig = ws->var("nsig");
RooRealVar* nBbkg = ws->var("nBbkg");
RooRealVar* nbkg = ws->var("nbkg");
RooRealVar* nbkg2 = ws->var("nbkg2");
RooDataSet* data = (RooDataSet*) ws->data("data");
// fit the model to the data.
model->fitTo(*data, Extended() );
RooMsgService::instance().setSilentMode(true);
// Now we use the SPlot class to add SWeights to our data set
// based on our model and our yield variables
RooStats::SPlot* sData = new RooStats::SPlot("sData","An SPlot",
*data, model, RooArgList(*nsig,*nBbkg,*nbkg,*nbkg2) );
// Check that our weights have the desired properties
std::cout << "Check SWeights:" << std::endl;
std::cout << std::endl << "Yield of sig is "
<< nsig->getVal() << ". From sWeights it is "
<< sData->GetYieldFromSWeight("nsig") << std::endl;
std::cout << std::endl << "Yield of Bbkg is "
<< nBbkg->getVal() << ". From sWeights it is "
<< sData->GetYieldFromSWeight("nBbkg") << std::endl;
std::cout << std::endl << "Yield of bkg is "
<< nbkg->getVal() << ". From sWeights it is "
<< sData->GetYieldFromSWeight("nbkg") << std::endl;
std::cout << std::endl << "Yield of bkg2 is "
<< nbkg2->getVal() << ". From sWeights it is "
<< sData->GetYieldFromSWeight("nbkg2") << std::endl;
cout << endl; cout << endl; cout << endl;
float sum20=0;
float sum50=0;
float sum100=0;
float sum200=0;
float sum300=0;
float sum600=0;
float sum900=0;
float sum1200=0;
float total=0;
// saving weights into a file
ofstream myfile;
myfile.open ("weights.txt");
// plot the weight event by event with the Sum of events values as cross-check
for(Int_t i=0; i < data->numEntries(); i++) {
//myfile << sData->GetSWeight(i,"nsig") << " " << sData->GetSWeight(i,"nBbkg") << " " << sData->GetSWeight(i,"nbkg") << " " << sData->GetSWeight(i,"nbkg2") << endl;
//myfile << sData->GetSWeight(i,"nsig") <<endl;
myfile << (unsigned int) data->get(i)->getRealValue("run")
<< " " << (unsigned int) data->get(i)->getRealValue("event")
<< " " << (float) data->get(i)->getRealValue("FourMu_Mass")
<< " " << sData->GetSWeight(i,"nsig")
<< endl;
// std::cout << "nsig Weight " << sData->GetSWeight(i,"nsig")
// << " nBbkg Weight " << sData->GetSWeight(i,"nBbkg")
// << " nbkg Weight " << sData->GetSWeight(i,"nbkg")
// << " nbkg2 Weight " << sData->GetSWeight(i,"nbkg2")
// << " Total Weight " << sData->GetSumOfEventSWeight(i)
// << std::endl;
total+=sData->GetSWeight(i,"nsig");
if(i<20) sum20+=sData->GetSWeight(i,"nsig");
if(i<50) sum50+=sData->GetSWeight(i,"nsig");
if(i<100) sum100+=sData->GetSWeight(i,"nsig");
if(i<200) sum200+=sData->GetSWeight(i,"nsig");
if(i<300) sum300+=sData->GetSWeight(i,"nsig");
if(i<600) sum600+=sData->GetSWeight(i,"nsig");
if(i<900) sum900+=sData->GetSWeight(i,"nsig");
if(i<1200) sum1200+=sData->GetSWeight(i,"nsig");
}
myfile.close();
std::cout << std::endl;
std::cout<<"Sum of the sWeights is: "<<total<<std::endl;
std::cout<<"Sum of the first 20 sWeights is: "<<sum20<<std::endl;
std::cout<<"Sum of the first 50 sWeights is: "<<sum50<<std::endl;
std::cout<<"Sum of the first 100 sWeights is: "<<sum100<<std::endl;
std::cout<<"Sum of the first 200 sWeights is: "<<sum200<<std::endl;
std::cout<<"Sum of the first 300 sWeights is: "<<sum300<<std::endl;
std::cout<<"Sum of the first 600 sWeights is: "<<sum600<<std::endl;
std::cout<<"Sum of the first 900 sWeights is: "<<sum900<<std::endl;
std::cout<<"Sum of the first 1200 sWeights is: "<<sum1200<<std::endl;
std::cout<<"Total # of events: "<<data->numEntries()<<std::endl;
// import this new dataset with sWeights
std::cout << "import new dataset with sWeights" << std::endl;
ws->import(*data, Rename("dataWithSWeights"));
//.........这里部分代码省略.........
示例12: drawHistoWithCurve
void drawHistoWithCurve( DrawBase* db, const std::string& data_dataset, const std::string& PUType, const std::string& data_mc, const std::string& flags, int nbtags, const std::string& leptType, std::string suffix ) {
if( suffix!="" ) suffix = "_" + suffix;
TH1F::AddDirectory(kTRUE);
// get histograms:
std::vector< TH1D* > lastHistos_data = db->get_lastHistos_data();
std::vector< TH1D* > lastHistos_mc = db->get_lastHistos_mc();
std::vector< TH1D* > lastHistos_mc_superimp = db->get_lastHistos_mc_superimp();
TH1D* h1_data = new TH1D(*(lastHistos_data[0]));
float xMin = (db->get_xAxisMin()!=9999.) ? db->get_xAxisMin() : h1_data->GetXaxis()->GetXmin();
float xMax = (db->get_xAxisMax()!=9999.) ? db->get_xAxisMax() : h1_data->GetXaxis()->GetXmax();
// create data graph (poisson asymm errors):
TGraphAsymmErrors* graph_data_poisson = new TGraphAsymmErrors(0);
graph_data_poisson = fitTools::getGraphPoissonErrors(h1_data);
graph_data_poisson->SetMarkerStyle(20);
THStack* mc_stack = new THStack();
for( unsigned ihisto=0; ihisto<lastHistos_mc.size(); ++ihisto )
mc_stack->Add(lastHistos_mc[lastHistos_mc.size()-ihisto-1]);
std::string sherpa_suffix = (use_sherpa) ? "_sherpa" : "";
// open fit results file:
char fitResultsFileName[200];
sprintf( fitResultsFileName, "fitResultsFile_%s_%dbtag_ALL_PU%s_fit%s%s%s.root", data_dataset.c_str(), nbtags, PUType.c_str(), data_mc.c_str(), sherpa_suffix.c_str(), flags.c_str() );
TFile* fitResultsFile = TFile::Open(fitResultsFileName);
// get bg workspace:
char workspaceName[200];
sprintf( workspaceName, "fitWorkspace_%dbtag", nbtags );
RooWorkspace* bgws = (RooWorkspace*)fitResultsFile->Get(workspaceName);
char fitResultsName[200];
sprintf( fitResultsName, "fitResults_%dbtag_decorr", nbtags );
RooFitResult* bgfr = (RooFitResult*)fitResultsFile->Get(fitResultsName);
// get mZZ variable:
RooRealVar* CMS_hzz2l2q_mZZ = (RooRealVar*)bgws->var("CMS_hzz2l2q_mZZ");
// get bg shape:
RooAbsPdf* background = (RooAbsPdf*)bgws->pdf("background_decorr");
std::string leptType_cut;
if( leptType=="MU" ) leptType_cut = "&& leptType==0";
if( leptType=="ELE" ) leptType_cut = "&& leptType==1";
SidebandFitter* sf = new SidebandFitter( data_dataset, PUType, data_mc, flags );
// get bg normalization:
//float expBkg = sf->get_backgroundNormalization( nbtags, leptType, "DATA" );
std::pair<Double_t, Double_t> bgNormAndError = sf->get_backgroundNormalizationAndError( nbtags, leptType, "DATA" );
float expBkg = bgNormAndError.first;
float expBkg_err = bgNormAndError.second;
//TTree* treeSidebandsDATA_alphaCorr = (TTree*)fitResultsFile->Get("sidebandsDATA_alpha");
//TH1D* h1_mZZ_sidebands_alpha = new TH1D("mZZ_sidebands_alpha", "", 65, xMin, xMax);
//char sidebandsCut_alpha[500];
//sprintf(sidebandsCut_alpha, "eventWeight_alpha*(isSidebands && nBTags==%d %s)", nbtags, leptType_cut.c_str());
//treeSidebandsDATA_alphaCorr->Project("mZZ_sidebands_alpha", "CMS_hzz2l2q_mZZ", sidebandsCut_alpha);
//float expBkg = h1_mZZ_sidebands_alpha->Integral();
RooPlot *plot_MCbkg = CMS_hzz2l2q_mZZ->frame(xMin,xMax,(int)(xMax-xMin)/h1_data->GetXaxis()->GetBinWidth(1));
background->plotOn(plot_MCbkg,RooFit::Normalization(expBkg));
if( suffix == "_turnOnZOOM" && !QUICK_ ) {
background->plotOn(plot_MCbkg,RooFit::VisualizeError(*bgfr,2.0,kFALSE),RooFit::FillColor(kYellow), RooFit::Normalization(expBkg));
background->plotOn(plot_MCbkg,RooFit::VisualizeError(*bgfr,1.0,kFALSE),RooFit::FillColor(kGreen), RooFit::Normalization(expBkg));
background->plotOn(plot_MCbkg,RooFit::Normalization(expBkg));
background->plotOn(plot_MCbkg,RooFit::LineStyle(2),RooFit::Normalization(expBkg+expBkg_err));
background->plotOn(plot_MCbkg,RooFit::LineStyle(2),RooFit::Normalization(expBkg-expBkg_err));
}
TF1* f1_bgForLegend = new TF1("bgForLegend", "[0]");
f1_bgForLegend->SetLineColor(kBlue);
f1_bgForLegend->SetLineWidth(3);
TH2D* h2_axes = new TH2D("axes", "", 10, xMin, xMax, 10, 0., 1.3*h1_data->GetMaximum());
char yTitle[200];
sprintf( yTitle, "Events / (%.0f GeV)", h1_data->GetXaxis()->GetBinWidth(1) );
h2_axes->SetYTitle(yTitle);
if( leptType=="MU" )
h2_axes->SetXTitle("m_{#mu#mujj} [GeV]");
else if( leptType=="ELE" )
h2_axes->SetXTitle("m_{eejj} [GeV]");
else
h2_axes->SetXTitle("m_{ZZ} [GeV]");
float legend_xMin = 0.38;
//.........这里部分代码省略.........
示例13: StandardHypoTestDemo
//.........这里部分代码省略.........
ropl = new RatioOfProfiledLikelihoodsTestStat(*bModel->GetPdf(), *sbModel->GetPdf(), sbModel->GetSnapshot());
ropl->SetSubtractMLE(false);
if (testStatType == 3) profll->SetOneSidedDiscovery(1);
profll->SetPrintLevel(printLevel);
// profll.SetReuseNLL(mOptimize);
// slrts.SetReuseNLL(mOptimize);
// ropl.SetReuseNLL(mOptimize);
AsymptoticCalculator::SetPrintLevel(printLevel);
HypoTestCalculatorGeneric * hypoCalc = 0;
// note here Null is B and Alt is S+B
if (calcType == 0) hypoCalc = new FrequentistCalculator(*data, *sbModel, *bModel);
else if (calcType == 1) hypoCalc= new HybridCalculator(*data, *sbModel, *bModel);
else if (calcType == 2) hypoCalc= new AsymptoticCalculator(*data, *sbModel, *bModel);
if (calcType == 0)
((FrequentistCalculator*)hypoCalc)->SetToys(ntoys, ntoys/nToysRatio);
if (calcType == 1)
((HybridCalculator*)hypoCalc)->SetToys(ntoys, ntoys/nToysRatio);
if (calcType == 2 ) {
if (testStatType == 3) ((AsymptoticCalculator*) hypoCalc)->SetOneSidedDiscovery(true);
if (testStatType != 2 && testStatType != 3)
Warning("StandardHypoTestDemo","Only the PL test statistic can be used with AsymptoticCalculator - use by default a two-sided PL");
}
// check for nuisance prior pdf in case of nuisance parameters
if (calcType == 1 && (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters() )) {
RooAbsPdf * nuisPdf = 0;
if (nuisPriorName) nuisPdf = w->pdf(nuisPriorName);
// use prior defined first in bModel (then in SbModel)
if (!nuisPdf) {
Info("StandardHypoTestDemo","No nuisance pdf given for the HybridCalculator - try to deduce pdf from the model");
if (bModel->GetPdf() && bModel->GetObservables() )
nuisPdf = RooStats::MakeNuisancePdf(*bModel,"nuisancePdf_bmodel");
else
nuisPdf = RooStats::MakeNuisancePdf(*sbModel,"nuisancePdf_sbmodel");
}
if (!nuisPdf ) {
if (bModel->GetPriorPdf()) {
nuisPdf = bModel->GetPriorPdf();
Info("StandardHypoTestDemo","No nuisance pdf given - try to use %s that is defined as a prior pdf in the B model",nuisPdf->GetName());
}
else {
Error("StandardHypoTestDemo","Cannnot run Hybrid calculator because no prior on the nuisance parameter is specified or can be derived");
return;
}
}
assert(nuisPdf);
Info("StandardHypoTestDemo","Using as nuisance Pdf ... " );
nuisPdf->Print();
const RooArgSet * nuisParams = (bModel->GetNuisanceParameters() ) ? bModel->GetNuisanceParameters() : sbModel->GetNuisanceParameters();
RooArgSet * np = nuisPdf->getObservables(*nuisParams);
if (np->getSize() == 0) {
Warning("StandardHypoTestDemo","Prior nuisance does not depend on nuisance parameters. They will be smeared in their full range");
}
delete np;
((HybridCalculator*)hypoCalc)->ForcePriorNuisanceAlt(*nuisPdf);
((HybridCalculator*)hypoCalc)->ForcePriorNuisanceNull(*nuisPdf);示例14: fitpeaks
//.........这里部分代码省略.........
RooCBShape *gauss1S2 = new RooCBShape ("gauss1S2", "FSR cb 1s",
*mass,*mean1S,*reso1S,*alpha,*npow);
//RooAddPdf *sig1S = new RooAddPdf ("sig1S","1S mass pdf",
// RooArgList(*gauss1S1,*gauss1S2),*sigmaFraction);
//mean->setVal(9.46);
//mean->setConstant(kTRUE);
sigma1->setVal(0);
sigma1->setConstant(kTRUE);
if (!fitMB) {
sigma2->setVal(width_); //fix the resolution
sigma2->setConstant(kTRUE);
}
/// Upsilon 2S
RooCBShape *gauss2S1 = new RooCBShape ("gauss2S1", "FSR cb 2s",
*mass,*mean2S,*sigma1,*alpha,*npow);
RooCBShape *gauss2S2 = new RooCBShape ("gauss2S2", "FSR cb 2s",
*mass,*mean2S,*reso2S,*alpha,*npow);
RooAddPdf *sig2S = new RooAddPdf ("sig2S","2S mass pdf",
RooArgList(*gauss2S1,*gauss2S2),*sigmaFraction);
/// Upsilon 3S
RooCBShape *gauss3S1 = new RooCBShape ("gauss3S1", "FSR cb 3s",
*mass,*mean3S,*sigma1,*alpha,*npow);
RooCBShape *gauss3S2 = new RooCBShape ("gauss3S2", "FSR cb 3s",
*mass,*mean3S,*reso3S,*alpha,*npow);
RooAddPdf *sig3S = new RooAddPdf ("sig3S","3S mass pdf",
RooArgList(*gauss3S1,*gauss3S2),*sigmaFraction);
/// Background
RooRealVar *bkg_a1 = new RooRealVar("bkg_{a1}", "background a1", 0, -2, 2);
RooRealVar *bkg_a2 = new RooRealVar("bkg_{a2}", "background a2", 0, -1, 1);
//RooRealVar *bkg_a3 = new RooRealVar("bkg_{a3}", "background a3", 0, -1, 1);
RooAbsPdf *bkgPdf = new RooChebychev("bkg","background",
*mass, RooArgList(*bkg_a1,*bkg_a2));
//bkg_a1->setVal(0);
//bkg_a1->setConstant(kTRUE);
//bkg_a2->setVal(0);
//bkg_a2->setConstant(kTRUE); //set constant for liner background
// only sideband region pdf, using RooPolynomial instead of RooChebychev for multiple ranges fit
RooRealVar *SB_bkg_a1 = new RooRealVar("SB bkg_{a1}", "background a1", 0, -1, 1);
RooRealVar *SB_bkg_a2 = new RooRealVar("SB bkg_{a2}", "background a2", 0, -1, 1);
RooAbsPdf *SB_bkgPdf = new RooPolynomial("SB_bkg","side-band background",
*mass, RooArgList(*SB_bkg_a1,*SB_bkg_a2));
//SB_bkg_a1->setVal(0);
//SB_bkg_a1->setConstant(kTRUE);
//SB_bkg_a2->setVal(0);
//SB_bkg_a2->setConstant(kTRUE);
/// Combined pdf
int nt = 100000;
//bool fitfraction = true;
RooRealVar *nbkgd = new RooRealVar("N_{bkg}","nbkgd",nt*0.75,0,10*nt);
RooRealVar *SB_nbkgd = new RooRealVar("SB N_{bkg}","SB_nbkgd",nt*0.75,0,10*nt);
RooRealVar *nsig1f = new RooRealVar("N_{#Upsilon(1S)}","nsig1S",nt*0.25,0,10*nt);
/*
//use the YIELDs of 2S and 3S as free parameters
RooRealVar *nsig2f = new RooRealVar("N_{#Upsilon(2S)}","nsig2S", nt*0.25,-1*nt,10*nt);
RooRealVar *nsig3f = new RooRealVar("N_{#Upsilon(3S)}","nsig3S", nt*0.25,-1*nt,10*nt);
*/
//use the RATIOs of 2S and 3S as free parameters
RooRealVar *f2Svs1S = new RooRealVar("N_{2S}/N_{1S}","f2Svs1S",0.21,-0.1,1);
//RooRealVar *f3Svs1S = new RooRealVar("N_{3S}/N_{1S}","f3Svs1S",0.0,-0.1,0.5);
RooRealVar *f23vs1S = new RooRealVar("N_{2S+3S}/N_{1S}","f23vs1S",0.45,-0.1,1);
RooFormulaVar *nsig2f = new RooFormulaVar("nsig2S","@0*@1", RooArgList(*nsig1f,*f2Svs1S));示例15: buildAndFitModels
void buildAndFitModels(TDirectory *fout, RooWorkspace *wspace, RooRealVar &x, std::string proc="Zvv"){
// Build and fit the model for the Zvv/Wlv background
RooAbsPdf *pdfZvv = wspace->pdf(doubleexp(wspace,x,Form("%s_control",proc.c_str())));
RooAbsPdf *pdfZvv_mc = wspace->pdf(doubleexp(wspace,x,Form("%s_control_mc",proc.c_str())));
RooAbsPdf *pdfZvv_background_mc = wspace->pdf(doubleexp(wspace,x,Form("%s_control_bkg_mc",proc.c_str())));
pdfZvv_mc->Print("v");
// Fit control region MC
std::cout << " Fit for control MC " << Form("%s_control_mc",proc.c_str())<< std::endl;
RooFitResult *fit_res_control_mc = pdfZvv_mc->fitTo(*(wspace->data(Form("%s_control_mc",proc.c_str()))),RooFit::Save(1),RooFit::SumW2Error(false));
fout->cd(); fit_res_control_mc->SetName(Form("fitResult_%s_control_mc",proc.c_str())); fit_res_control_mc->Write();
std::cout << " Fit for background MC " << Form("%s_control_bkg_mc",proc.c_str()) << std::endl;
// Fit background MC and then fix it
pdfZvv_background_mc->fitTo(*(wspace->data(Form("%s_control_bkg_mc",proc.c_str()))),RooFit::SumW2Error(true));
freezeParameters(pdfZvv_background_mc->getParameters(RooArgSet(x)));
// Now fit the Zvv Data
//RooRealVar frac_contamination_Zvv(Form("frac_contamination_%s",proc.c_str()),Form("frac_contamination_%s",proc.c_str()),0,1);
double nbkgcont = wspace->data(Form("%s_control_bkg_mc",proc.c_str()))->sumEntries();
double ncont = wspace->data(Form("%s_control",proc.c_str()))->sumEntries()-nbkgcont;
RooRealVar num_contamination_Zvv(Form("num_contamination_%s",proc.c_str()),Form("num_contamination_%s",proc.c_str()),nbkgcont,0,10E10);
num_contamination_Zvv.setConstant();
RooRealVar num_Zvv(Form("num_%s",proc.c_str()),Form("num_%s",proc.c_str()),ncont,0,10E10);
num_Zvv.setConstant(true);// freeze the n_data now
RooAddPdf modelZvv(Form("model_%s_control",proc.c_str()),Form("model_%s_control",proc.c_str()),RooArgList(*pdfZvv_background_mc,*pdfZvv),RooArgList(num_contamination_Zvv,num_Zvv));
std::cout << " Fit for control Data " << Form("%s_control",proc.c_str()) << std::endl;
RooFitResult *fit_res_control = modelZvv.fitTo(*(wspace->data(Form("%s_control",proc.c_str()))),RooFit::Save(1));
fout->cd(); fit_res_control->SetName(Form("fitResult_%s_control",proc.c_str())); fit_res_control->Write();
// Find the scale of ndata/nmc to normalize the yields
double nmontecarlo = wspace->data(Form("%s_control_mc",proc.c_str()))->sumEntries();
double ndata = num_Zvv.getVal();
RooRealVar scalef(Form("scalef_%s",proc.c_str()),"scalef",ndata/nmontecarlo);
// uncomment make this ONLY a shape correction!
// scalef.setVal(1);
scalef.setConstant(true);
std::cout << proc.c_str() << " N Control Data == " << ndata << std::endl;
std::cout << proc.c_str() << " N Control MC == " << nmontecarlo << std::endl;
// Plot the fits
TCanvas can_datafit(Form("%s_datafit",proc.c_str()),"Data Fit",800,600);
RooPlot *pl = x.frame();
(wspace->data(Form("%s_control_bkg_mc",proc.c_str())))->plotOn(pl,RooFit::MarkerStyle(kOpenCircle));
(wspace->data(Form("%s_control",proc.c_str())))->plotOn(pl);
modelZvv.plotOn(pl);
modelZvv.paramOn(pl);
//pdfZvv_background_mc->plotOn(pl,RooFit::LineStyle(2));
pl->Draw();
fout->cd();can_datafit.Write();
TCanvas can_mcfit(Form("%s_mcfit",proc.c_str()),"MC Fit",800,600);
RooPlot *plmc = x.frame();
(wspace->data(Form("%s_control_mc",proc.c_str())))->plotOn(plmc,RooFit::MarkerColor(kBlack));
pdfZvv_mc->plotOn(plmc,RooFit::LineStyle(1),RooFit::LineColor(2));
pdfZvv_mc->paramOn(plmc);
plmc->Draw();
fout->cd();can_mcfit.Write();
TCanvas can_mcdatafit(Form("%s_mcdatafit",proc.c_str()),"MC and Data Fits",800,600);
RooPlot *plmcdata = x.frame();
pdfZvv_mc->plotOn(plmcdata,RooFit::LineColor(2),RooFit::Normalization(nmontecarlo));
pdfZvv->plotOn(plmcdata,RooFit::Normalization(ndata));
plmcdata->Draw();
fout->cd();can_mcdatafit.Write();
// Import the correction and the models
RooFormulaVar ratio(Form("ratio_%s",proc.c_str()),"@0*@1/@2",RooArgList(scalef,*pdfZvv,*pdfZvv_mc));
wspace->import(ratio);
wspace->import(num_Zvv);
TCanvas can_ra(Form("%s_ratio",proc.c_str()),"MC Fit",800,600);
RooPlot *plra = x.frame();
ratio.plotOn(plra,RooFit::LineStyle(1));
plra->Draw();
fout->cd();can_ra.Write();
}