Analisi Statistica dei Dati per HEP (Laboratorio)

1. Analisi Statistica dei Dati per HEP (Laboratorio) • 2013-03-25 • Elementi di C++ • Introduzione a ROOT • 2013-04-11, Laboratorio Informatico • ROOT warm up • 2013-04-17, Laboratorio Informatico • Introduzione a RooFit • Primo esercizio con RooFit • 2013-05-02, Laboratorio Informatico • RooFit (Workspace, Factory, Composite Model) • 2013-05-16, Laboratorio Informatico • Introduzione a RooStats • 2013-05-23, Laboratorio Informatico • RooStats , TMVA

2. RECAP: Esercizio RooStats [1] roostats_ex1.C • Si modifichi il modello:Specificareicomponenti del modello per i tool statistici di roostat: osservabile e parametro di interesse. Utilizzareilnumero di eventi di segnale come unicoparametro di interesse. Fissarecostantituttiglialtriparametri del modello.Importare la configurazionenel workspace e salvaresu file. • Si modifichil’uso del modello:- leggereilmodelConfig dal workspace esempio: ModelConfig* mc = (ModelConfig*) w.obj("ModelConfig");- calcolare un Confidence Interval utilizzandoilProfileLikelihoodCalculator- Disegnareilprofilodella likelihood e sovrapporrel’intervallo- calcolare la discovery significance utilizzandoilprofilelikelihoodcalculator come test di ipotesi- scriveresulla console ilimitidell’intervallo e la significatività [2] roostats_ex1.C • Aggiungere l’intervallo calcolato con Feldman-Cousinssuggerimento : guardare il codice in :$ROOTSYS/tuturials/roostats/IntervalExamples.C

3. voidmakeModel(RooWorkspace& w ) { // Construct model here w.factory("Exponential:bkg(x[110,120], expr('(-1./tau)',tau[10,-1000,1000]) )"); w.factory("Gaussian:sig(x, mean[115,105,120], sigma[1.0,0,10])"); w.factory("SUM:model( Ns[100,0,120]*sig, Nb[1000,0,10000]*bkg)"); // Getting variables and functions out of a workspace RooAbsPdf * model = w.pdf("model"); RooRealVar * x = w.var("x"); // the observable RooRealVar * Ns = w.var("Ns"); // the parameter x->setBins(30); // set constant parameters w.var("mean")->setConstant(true); w.var("sigma")->setConstant(true); w.var("tau")->setConstant(true); w.var("Nb")->setConstant(true); //specify components of model for statistical tools ModelConfigmodelConfig("ModelConfig"); modelConfig.SetWorkspace(w); //set components using the name of ws objects modelConfig.SetPdf( "model" ); modelConfig.SetObservables("x"); modelConfig.SetParametersOfInterest("Ns"); //import modelConfig into workspace too w.import(modelConfig); // generate the data (nsig = 50, nbkg=1000) w.var("Ns")->setVal(50); w.var("Nb")->setVal(1000); // use fixed random numbers for reproducibility RooRandom::randomGenerator()->SetSeed(111); intN = 1100 ; // will generate accordint to total S+B events RooDataSet* data = model->generate( *x ); data->Print(); // add data to the workspace data->SetName("data"); w.import(*data); // write the workspace in the file constchar * fileName = "model.root"; w.writeToFile(fileName,true); } non metto N

4. voiduseModel(RooWorkspace& w ) { // Getting variables and functions out of a workspace RooAbsPdf * model = w.pdf("model"); RooRealVar * x = w.var("x"); // the observable RooRealVar * Ns = w.var("Ns"); // the parameter RooDataSet * data = w.data("data"); // the data // get the modelConfig out of the ModelConfig* mc = (ModelConfig*) w.obj("ModelConfig"); // create the class using data and model ProfileLikelihoodCalculatorplc(*data, *mc); // set the confidencelevel plc.SetConfidenceLevel(0.683); // compute the interval LikelihoodInterval* interval = plc.GetInterval(); double lowerLimit = interval->LowerLimit(*Ns); double upperLimit = interval->UpperLimit(*Ns); // plot the interval LikelihoodIntervalPlot plot(interval); plot.Draw(); // Calculate the significance: SL2 estimator Ns->setVal(0); plc.SetNullParameters(*Ns); HypoTestResult* hypotest = plc.GetHypoTest(); double SL2 = hypotest->Significance(); cout << "Limits: [" << lowerLimit << " - " << upperLimit << "]" << endl; cout << "Significance (SL2): " << SL2 << endl; } Chiama il ProfileLikelihoodCalculator Passando gli argomenti (dati, modelconfig)

5. 68% Limits: [27.9087 - 67.6207] Significance (SL2): 2.50742 68% C.L. Ns RicordarecheNbera fissato costante, così come gli altri parametri del modello (tau, mass, sigma).

6. voiduseModel(RooWorkspace& w ) { // Getting variables and functions out of a workspace RooAbsPdf * model = w.pdf("model"); RooRealVar * x = w.var("x"); // the observable RooRealVar * Ns = w.var("Ns"); // the parameter RooDataSet * data = w.data("data"); // the data // get the modelConfig out of the ModelConfig* mc = (ModelConfig*) w.obj("ModelConfig"); FeldmanCousinsfc(*data, *modelConfig); fc.SetConfidenceLevel( 0.95 ); fc.SetNBins(40); // number of points to test per parameter fc.UseAdaptiveSampling(true); // make it go faster // The PDF could be extended and this could be removed // fc.FluctuateNumDataEntries(false); // Proof // ProofConfig pc(*wspace, 4, "workers=4", kFALSE); // proof-lite //ProofConfig pc(w, 8, "localhost"); // proof cluster at "localhost" // ToyMCSampler* toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler(); // toymcsampler->SetProofConfig(&pc); // enableproof PointSetInterval* interval = (PointSetInterval*) fc.GetInterval(); std::cout << “Feldman-Cousins Limits: ["<< interval->LowerLimit(*Ns) << " , " << interval->UpperLimit(*Ns) << "]" << endl; } ESERCIZIO 2 Intervallo Feldman-Cousins Chiama il FeldmanCousinscalculator Passando sempre gli argomenti (dati, modelconfig)

7. === Using the following for ModelConfig === Observables: RooArgSet:: = (x) Parameters of Interest: RooArgSet:: = (Ns) PDF: RooAddPdf::model[ Ns * sig + Nb * bkg ] FeldmanCousins: ntoys per point: adaptive FeldmanCousins: nEvents per toywillfluctuateaboutexpectation FeldmanCousins: Model has no nuisanceparameters FeldmanCousins: # points to test = 40 FeldmanCousinsIntervalis[16.5 , 85.5 ]ProfileLikeLihoodIntervalis[9.94475 , 87.7451] Significance(SL2):2.50742 95% C.L. Ns 95% C.L. RicordarecheNbera fissato costante, così come gli altri parametri del modello (tau, mass, sigma).

8. voidmakeModel(RooWorkspace& w ) { // Construct model here w.factory("Exponential:bkg(x[110,120], expr('(-1./tau)',tau[10,-1000,1000]) )"); w.factory("Gaussian:sig(x, mean[115,105,120], sigma[1.0,0,10])"); w.factory("SUM:model( Ns[100,0,120]*sig, Nb[1000,0,10000]*bkg)"); // Getting variables and functions out of a workspace RooAbsPdf * model = w.pdf("model"); RooRealVar * x = w.var("x"); // the observable RooRealVar * Ns = w.var("Ns"); // the parameter x->setBins(30); // set constant parameters w.var("mean")->setConstant(true); w.var("sigma")->setConstant(true); // define set of nuisanceparameters w.defineSet("nuisParams","tau,Nb"); //specify components of model for statistical tools ModelConfigmodelConfig("ModelConfig"); modelConfig.SetWorkspace(w); //set components using the name of ws objects modelConfig.SetPdf( "model" ); modelConfig.SetObservables("x"); modelConfig.SetParametersOfInterest("Ns"); modelConfig.SetNuisanceParameters(*w.set("nuisParams")); //import modelConfig into workspace too w.import(modelConfig); VARIAZIONEIntroduciamo NuisanceParameters

9. FeldmanCousinsIntervalis[16.5 , 100.5 ] ProfileLikeLihoodIntervalis [10.1149 - 106.863] Significance (SL2): 2.3894 95% C.L. Ns 95% C.L.

10. VARIAZIONEIntervallo Bayesian // example use of BayesianCalculator // now we also need to specify a prior in the ModelConfig wspace->factory("Uniform::prior(mu)"); modelConfig->SetPriorPdf(*wspace->pdf("prior")); // exampleusage of BayesianCalculator BayesianCalculatorbc(*data, *modelConfig); bc.SetConfidenceLevel( confidenceLevel); SimpleInterval* bcInt = bc.GetInterval(); Chiama il BayesianCalculatorPassando sempre gli argomenti (dati, modelconfig)

11. We create the class and configure it. Since the numerical integration can be tricky, it is better in this case to set a reasonable range from the nuisance parameters which will be integrated to obtain the marginalised posterior function. A sensible choice is to use for example an interval using 10-sigma around the best fit values for the nuisance parameters, RooRealVar *nuisPar = w->var("..."); nuisPar->setRange(nuisPar->getVal() - 10*nuisPar->getError(), nuisPar->getVal() + 10* nuisPar->getError() ) ; It is also recommended to use this class on a model which has a binned dataset or not having too many events (<~ 100), otherwise it will take a long time to obtain the result. We create the class can set the type of interval (central, upper/lower limit or shortest) and the type of integration method BayesianCalculatorbayesianCalc(*data,*mc); bayesianCalc.SetConfidenceLevel(0.683); // 68% interval // set the type of interval (not really needed for central which is the default) bayesianCalc.SetLeftSideTailFraction(0.5); // for central interval //bayesianCalc.SetLeftSideTailFraction(0.); // for upper limit //bayesianCalc.SetShortestInterval(); // for shortest interval // set the integration type (not really needed for the default ADAPTIVE) // possible alternative values are "VEGAS" , "MISER", or "PLAIN" (MC integration from libMathMore) // "TOYMC" (toy MC integration, work when nuisances exist and they have a constraints pdf) TStringintegrationType = ""; // this is needed if using TOYMC if (integrationType.Contains("TOYMC") ) { RooAbsPdf * nuisPdf = RooStats::MakeNuisancePdf(*mc, "nuisance_pdf"); if (nuisPdf) bayesianCalc.ForceNuisancePdf(*nuisPdf); } bayesianCalc.SetIntegrationType(integrationType);

12. // compute interval by scanning the posterior function // it is done by default when computing shortest intervals bayesianCalc.SetScanOfPosterior(100); RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first(); SimpleInterval* interval = bayesianCalc.GetInterval(); if (!interval) { cout << "Error computing Bayesian interval - exit " << endl; return; } double lowerLimit = interval->LowerLimit(); double upperLimit = interval->UpperLimit(); cout << "\n68% interval on " <<firstPOI->GetName()<<" is : ["<< lowerLimit << ", "<< upperLimit <<"] "<<endl; // draw plot of posterior function RooPlot * plot = bayesianCalc.GetPosteriorPlot(); if (plot) plot->Draw();

13. TMVA