00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00014
00015
00016
00017
00018
00019
00020 #include "RooRarFit/rarVersion.hh"
00021
00022 #include "Riostream.h"
00023 #include <sstream>
00024 #include <map>
00025 #include <set>
00026 #include <vector>
00027 using namespace std;
00028
00029 #include <libgen.h>
00030 #include "TFile.h"
00031 #include "TTree.h"
00032 #include "TObjString.h"
00033 #include "TStopwatch.h"
00034
00035 #include "RooFitCore/Roo1DTable.hh"
00036 #include "RooFitCore/RooAbsPdf.hh"
00037 #include "RooFitCore/RooAddPdf.hh"
00038 #include "RooFitCore/RooChi2Var.hh"
00039 #include "RooFitCore/RooHistError.hh"
00040 #include "RooFitCore/RooArgList.hh"
00041 #include "RooFitCore/RooDataSet.hh"
00042 #include "RooFitCore/RooFitResult.hh"
00043 #include "RooFitCore/RooMCStudy.hh"
00044 #include "RooFitCore/RooNLLVar.hh"
00045 #include "RooFitCore/RooPlot.hh"
00046 #include "RooFitCore/RooProdPdf.hh"
00047 #include "RooFitCore/RooRandom.hh"
00048 #include "RooFitCore/RooRealVar.hh"
00049 #include "RooFitCore/RooSimPdfBuilder.hh"
00050 #include "RooFitCore/RooStringVar.hh"
00051 #include "RooFitCore/RooSuperCategory.hh"
00052 #include "RooFitCore/RooGlobalFunc.hh"
00053 #include "RooFitCore/RooPullVar.hh"
00054
00055 using namespace RooFit;
00056
00057 #include "RooRarFit/rarMinuit.hh"
00058 #include "RooRarFit/rarMLPdf.hh"
00059 #include "RooRarFit/rarNLL.hh"
00060 #include "RooRarFit/rarSPlot.hh"
00061
00062 #include "RooRarFit/rarMLFitter.hh"
00063 #include "RooRarFit/rarToyList.hh"
00064
00065 ClassImp(rarMLFitter)
00066 ;
00067
00068 TString rarMLFitter::_physCatStr="";
00069 TString rarMLFitter::_splitCatStr="";
00070 RooArgSet rarMLFitter::_splitCatSet;
00071 RooArgSet rarMLFitter::_splitCatSet2;
00072
00076 rarMLFitter::rarMLFitter()
00077 : rarCompBase(),
00078 _simBuilder(0), _simConfig(0), _theGen(0), _protGenLevel(0),
00079 _protDataset(0), _theToyParamGen(0), _theSPdf(0), _theBPdf(0),
00080 _toyID(0), _toyNexp(0)
00081 {
00082 init();
00083 }
00084
00097 rarMLFitter::rarMLFitter(const char *configFile, const char *configSec,
00098 const char *configStr, rarDatasets *theDatasets,
00099 RooDataSet *theData, const char*name,const char*title)
00100 : rarCompBase(configFile, configSec, configStr,
00101 theDatasets, theData, name, title, kFALSE),
00102 _simBuilder(0), _simConfig(0), _theGen(0), _protGenLevel(0),
00103 _protDataset(0), _theToyParamGen(0), _theSPdf(0), _theBPdf(0),
00104 _toyID(0), _toyNexp(0)
00105 {
00106 init();
00107 }
00108
00109 rarMLFitter::~rarMLFitter()
00110 {
00111 }
00112
00118 void rarMLFitter::init()
00119 {
00120
00121 if (!getControlBit("SimFit"))
00122 _thePdf=(RooAbsPdf*)_subPdfs.at(0);
00123 else if (""!=readConfStr("category", "", getVarSec())) {
00124 TString catStr=readConfStr("category", "", getVarSec());
00125 _category=(RooAbsCategoryLValue *)_fullObs->find(catStr);
00126 if (!_category) {
00127 cout<<"Can not find category named "<<catStr<<endl;
00128 exit(-1);
00129 }
00130 _thePdf=new RooSimultaneous(Form("the_%s",GetName()),
00131 _pdfType+" "+GetTitle(),
00132 _subPdfs, *_category);
00133 } else {
00134
00135 _simBuilder=new RooSimPdfBuilder(_subPdfs);
00136 _simConfig=_simBuilder->createProtoBuildConfig();
00137 _simConfig->readFromFile(_configFile, 0, getVarSec());
00138 cout<<"simPdfBuilder configs in the config file:"<<endl;
00139 _simConfig->Print("v");
00140
00141 setSpecialStr(_simConfig);
00142 cout<<"simPdfBuilder configs after reading in special strings:"<<endl;
00143 _simConfig->Print("v");
00144
00145
00146 _simBuilder->addSpecializations(getSpecialSet());
00147 cout<<"special Set"<<endl;
00148 getSpecialSet().Print("v");
00149 if (!_theData) {
00150 cout<<" No dataset defined to build simPdf."
00151 <<" Please add/edit"<<endl
00152 <<"fitData=<datasetName>"<<endl
00153 <<" in master pdf section ["<<getMasterSec()<<"]"<<endl;
00154 exit(-1);
00155 }
00156 _thePdf=(RooAbsPdf*)_simBuilder->buildPdf(*_simConfig, _theData, 0, kTRUE);
00157 if (!_thePdf) {
00158 cout<<" simPdf not built!!"<<endl
00159 <<" Please check your simPdf config in section ["
00160 <<getVarSec()<<"]"<<endl;
00161 exit(-1);
00162 }
00163
00164
00165 setSimPdf((RooSimultaneous*)_thePdf);
00166 }
00167 {
00168 Int_t xPdfs=_xPdfList.GetSize();
00169 createPdfs("preToyRandGenerators",&_xPdfList,&_preToyRandGenerators,
00170 _runSec);
00171 Int_t xPdfsAfter=_xPdfList.GetSize();
00172 if (xPdfsAfter>xPdfs) {
00173 for (Int_t i=xPdfs; i<xPdfsAfter; i++) {
00174 rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
00175 thePdf->setControlBit("noPdfFit", "pdfFit");
00176 thePdf->setControlBit("noPdfPlot", "pdfPlot");
00177 }
00178
00179 _theToyParamGen=new RooProdPdf
00180 (Form("paraRand_%s", GetName()),"paraRand",_preToyRandGenerators);
00181 }
00182 }
00183
00184 rarStrParser postEmbdRandObsParser=readConfStr("postEmbdRandObs","",_runSec);
00185 TString embdObsGensStr="";
00186 while (postEmbdRandObsParser.nArgs()>0) {
00187
00188 TString dataSrcName=postEmbdRandObsParser[0];
00189 postEmbdRandObsParser.Remove();
00190 if (_embdObsRandSet.find(dataSrcName)) {
00191 cout<<dataSrcName<<" has been defined in postEmbdRandObs"<<endl;
00192 exit(-1);
00193 }
00194 RooStringVar *srcRandStr=new RooStringVar(dataSrcName,dataSrcName,"",8192);
00195 _embdObsRandSet.addOwned(*srcRandStr);
00196 while (postEmbdRandObsParser.nArgs()>0) {
00197 TString obsName=postEmbdRandObsParser[0];
00198 postEmbdRandObsParser.Remove();
00199 srcRandStr->setVal(Form("%s %s",srcRandStr->getVal(), obsName.Data()));
00200
00201 if (_thePdf->getDependents(*_fullObs)->find(obsName)) continue;
00202 embdObsGensStr+=" "+obsName;
00203 break;
00204 }
00205 }
00206 if (""!=embdObsGensStr) {
00207 setConfStr("embdObsGensStr", embdObsGensStr);
00208 Int_t xPdfs=_xPdfList.GetSize();
00209 createPdfs("embdObsGensStr",&_xPdfList,&_embdObsGens);
00210 Int_t xPdfsAfter=_xPdfList.GetSize();
00211 if (xPdfsAfter>xPdfs) {
00212 for (Int_t i=xPdfs; i<xPdfsAfter; i++) {
00213 rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
00214 thePdf->setControlBit("noPdfFit", "pdfFit");
00215 thePdf->setControlBit("noPdfPlot", "pdfPlot");
00216 }
00217 }
00218 }
00219
00220 if (_protDataVars.getSize()>0) {
00221 cout<<"protDataVars (including conditionalObs) defined:"<<endl;
00222 _protDataVars.Print("v");
00223 }
00224 if (_conditionalObs.getSize()>0) {
00225 cout<<"conditionalObs defined:"<<endl;
00226 _conditionalObs.Print("v");
00227 }
00228
00229 {
00230 RooArgSet allProtEVars(getArgSet("protDataEVars", kTRUE));
00231 allProtEVars.add(getArgSet(readConfStr("protDataEVars","",_runSec),kFALSE));
00232 if (allProtEVars.getSize() > 0) {
00233 cout<<"protDataEVars defined:"<<endl;
00234 allProtEVars.Print("v");
00235 } else {
00236 cout << "No protDataEVars defined." << endl;
00237 }
00238 RooArgList allProtEVarList(allProtEVars);
00239 for (Int_t i=0; i<allProtEVarList.getSize(); i++) {
00240 RooAbsArg *theEVar=_fullObs->find(allProtEVarList[i].GetName());
00241 if (theEVar) _protDataEVars.add(*theEVar);
00242 }
00243 }
00244
00245 cout<<"done init of rarMLFitter for "<<GetName()<<endl<<endl;
00246 }
00247
00253 void rarMLFitter::setSpecialStr(RooArgSet *simConfigSet)
00254 {
00255 for (Int_t i=0; i<_nComp; i++) {
00256 rarMLPdf *model=(rarMLPdf*)_pdfList.At(i);
00257 RooStringVar* modelSplitStr=(RooStringVar*)
00258 simConfigSet->find(Form("the_%s",model->GetName()));
00259 modelSplitStr->setVal(Form("%s %s", modelSplitStr->getVal(),
00260 (model->getSpecialStr()).Data()));
00261
00262 TString splitVar=modelSplitStr->getVal();
00263 Int_t lIdx, rIdx;
00264 while (((lIdx=splitVar.First("["))>-1)&&
00265 ((rIdx=splitVar.First("]"))>-1)) {
00266 splitVar.Replace(lIdx, rIdx-lIdx+1, "");
00267 }
00268 modelSplitStr->setVal(splitVar.Data());
00269 }
00270 }
00271
00278 RooArgSet rarMLFitter::getSpecialSet(TString setName)
00279 {
00280 RooArgSet specialSet("special ArgSet for splitting");
00281 for (Int_t i=0; i<_nComp; i++) {
00282 specialSet.add(((rarMLPdf*)_pdfList.At(i))->getSpecialSet(setName));
00283 }
00284 return specialSet;
00285 }
00286
00298 void rarMLFitter::getSplitCoeffValues(RooArgList coeffList, TString valType,
00299 ostream &o)
00300 {
00301 for (Int_t i=0; i<coeffList.getSize(); i++) {
00302 RooRealVar *theCoeffFrac=(RooRealVar*)(&coeffList[i]);
00303 Double_t val=theCoeffFrac->getVal();
00304 if ("Frac"==valType) {
00305 o<<theCoeffFrac->GetName()<<" = "<<val<<endl;
00306 } else if ("Asym"==valType) {
00307 Double_t asym=2*val-1;
00308 Double_t err=2*theCoeffFrac->getError();
00309 Double_t errHi=2*theCoeffFrac->getAsymErrorHi();
00310 Double_t errLo=2*theCoeffFrac->getAsymErrorLo();
00311 o<<"Asym wrt "<<theCoeffFrac->GetName()<<endl
00312 <<" 1-2*("<<theCoeffFrac->GetName()<<") = "
00313 <<-asym<<" +/- "<<err<<" (+"<<errLo<<", -"<<errHi<<")"<<endl;
00314
00315
00316 }
00317 }
00318 }
00319
00322 TString rarMLFitter::getPhysCat()
00323 {
00324 if (""!=_physCatStr) return _physCatStr;
00325 rarStrParser physModelsParser=readConfStr("physModels", "", getMasterSec());
00326 _physCatStr=physModelsParser[0];
00327 _physCatStr.ReplaceAll(":", "");
00328 RooAbsCategory *theCat=(RooAbsCategory*)(getCats()->find(_physCatStr));
00329 if (!theCat) _physCatStr="";
00330 return _physCatStr;
00331 }
00332
00335 TString rarMLFitter::getSplitCats()
00336 {
00337 if (""!=_splitCatStr) return _splitCatStr;
00338
00339 getSplitCatSet();
00340
00341 Int_t nCat=_splitCatSet.getSize();
00342 if (nCat<=0) {
00343 _splitCatStr=" ";
00344 return _splitCatStr;
00345 }
00346 RooArgList splitCatList(_splitCatSet);
00347 for(Int_t i=0; i<nCat; i++) {
00348 _splitCatStr+=splitCatList[i].GetName();
00349 _splitCatStr+=" ";
00350 }
00351 _splitCatStr.Remove(_splitCatStr.Length()-1);
00352
00353 return _splitCatStr;
00354 }
00355
00358 RooArgSet rarMLFitter::getSplitCatSet()
00359 {
00360 if (_splitCatSet.getSize()>0) return _splitCatSet;
00361
00362 TString splitCats=readConfStr("splitCats", "", getMasterSec());
00363 TString physCat=getPhysCat();
00364 if (""!=physCat) splitCats.ReplaceAll(physCat, "");
00365 rarStrParser catsStrParser=splitCats;
00366 Int_t nCat=catsStrParser.nArgs();
00367 for (Int_t i=0; i<nCat; i++) {
00368 TString catName=catsStrParser[i];
00369
00370 Int_t lIdx, rIdx;
00371 while (((lIdx=catName.First("("))>-1)&&
00372 ((rIdx=catName.First(")"))>-1)) {
00373 catName.Replace(lIdx, rIdx-lIdx+1, "");
00374 }
00375 if (""==catName) continue;
00376 RooAbsCategory *theCat=(RooAbsCategory*)(getCats()->find(catName));
00377 if (!theCat) {
00378 cout<<"Can not find cat "<<catName<<endl
00379 <<"Please make sure the cat you are using exists"<<endl;
00380 exit(-1);
00381 }
00382
00383 lIdx=catsStrParser[i].First("(");
00384 rIdx=catsStrParser[i].First(")");
00385 if (lIdx<0 && rIdx<0) {
00386 _splitCatSet.add(*theCat);
00387 continue;
00388 }
00389 if (lIdx>-1 && rIdx<0) {
00390 cout<<" Invalid specif. w/ config splitCats in sec. "<<getMasterSec()
00391 <<endl;
00392 exit(-1);
00393 }
00394 if (lIdx<0 && rIdx>-1) {
00395 cout<<" Invalid specif. w/ config splitCats in sec. "<<getMasterSec()
00396 <<endl;
00397 exit(-1);
00398 }
00399
00400 TString catTypesStr=catsStrParser[i](lIdx+1, rIdx-lIdx-1);
00401 catTypesStr.ReplaceAll(",", " ");
00402 rarStrParser catTypesStrParser=catTypesStr;
00403 if (catTypesStrParser.nArgs()<1) continue;
00404
00405 RooCategory *thePCat=new RooCategory(theCat->GetName(),theCat->GetTitle());
00406
00407 for (Int_t k=0; k<catTypesStrParser.nArgs(); k++) {
00408
00409 const RooCatType *theType=theCat->lookupType(catTypesStrParser[k]);
00410 if (!theType) {
00411 cout <<"The cat type "<<catTypesStrParser[k]
00412 <<" does not exist in cat "<<theCat->GetName()<<endl;
00413 exit(-1);
00414 }
00415 thePCat->defineType(theType->GetName(), theType->getVal());
00416 }
00417 _splitCatSet.add(*thePCat);
00418 }
00419
00420
00421
00422
00423 return _splitCatSet;
00424 }
00425
00430 RooAbsCategory *rarMLFitter::getSplitCat(RooArgSet &splitCatSet,
00431 TString catName)
00432 {
00433 TString retName="";
00434 RooAbsCategory *retVal(0);
00435 RooArgSet thisSet;
00436 rarStrParser catNameParser=catName;
00437 while (catNameParser.nArgs()>0) {
00438 TString theName=catNameParser[0];
00439 catNameParser.Remove();
00440 if (!splitCatSet.find(theName)) {
00441 cout<<" Can not find "<<theName<<" in split cats "<<endl;
00442 return retVal;
00443 }
00444 thisSet.add(*splitCatSet.find(theName));
00445 retName+=theName+" ";
00446 }
00447 retName.Remove(retName.Length()-1);
00448 retName.ReplaceAll(" ", "_");
00449
00450 splitCatSet.remove(thisSet);
00451
00452
00453 if (1==thisSet.getSize()) {
00454 retVal=(RooAbsCategory*)thisSet.find(retName);
00455 return retVal;
00456 }
00457
00458 retVal=(RooAbsCategory*)_splitCatSet2.find(retName);
00459 if (retVal) return retVal;
00460
00461 retVal=new RooSuperCategory(retName, retName, thisSet);
00462
00463 _splitCatSet2.add(*retVal);
00464
00465 return retVal;
00466 }
00467
00473 TString rarMLFitter::getRootFileName(TString aType, TString configToken)
00474 {
00475 TString rootFile;
00476 if (("yes"==configToken)||("default"==configToken))
00477 rootFile=getFullFileName(_resultDir,aType,_runSec,"none", getMasterSec());
00478 else {
00479 char *bName=basename((char*)configToken.Data());
00480 rootFile=_resultDir+"/"+bName;
00481 }
00482 if (!rootFile.EndsWith(".root")) rootFile+=".root";
00483 if (_toyID) {
00484 rootFile.Replace(rootFile.Length()-5, 5, ".%03d.root");
00485 rootFile=Form(rootFile.Data(), _toyID);
00486 }
00487
00488 return rootFile;
00489 }
00490
00497 TString rarMLFitter::getParamFileName(TString aType, TString configToken,
00498 TString dirName)
00499 {
00500 TString paramFile;
00501 TString name="none";
00502 TString dsName="";
00503 TString msName="";
00504 TString cfName=_configFile;
00505 rarStrParser paramIDParser=configToken;
00506 while (paramIDParser.nArgs()>1) {
00507 char flag=paramIDParser[0][0];
00508 paramIDParser.Remove();
00509 switch (flag) {
00510 case 'F' :
00511 cfName=paramIDParser[0];
00512 paramIDParser.Remove();
00513 break;
00514 case 'D' :
00515 dsName=paramIDParser[0];
00516 paramIDParser.Remove();
00517 break;
00518 case 'C' :
00519 msName=paramIDParser[0];
00520 paramIDParser.Remove();
00521 break;
00522 case 'A' :
00523 aType=paramIDParser[0];
00524 paramIDParser.Remove();
00525 break;
00526 case 'N' :
00527 name=paramIDParser[0];
00528 paramIDParser.Remove();
00529 break;
00530 }
00531 }
00532
00533 if (""==dirName) dirName=_paramDir;
00534 return getFullFileName(dirName, aType, name, dsName, msName, cfName);
00535 }
00536
00543 void rarMLFitter::paramFileIO(RooArgSet params, TString paramFile, Bool_t In)
00544 {
00545
00546 RooArgSet corrCoeffs(getCorrCoeffs());
00547 TString paramOrder=readConfStr("outParamOrder", "ascend", getMasterSec());
00548 if (("ascend"==paramOrder) || ("descend"==paramOrder)) {
00549 Bool_t reverse=("descend"==paramOrder);
00550 RooArgList paramList(corrCoeffs);
00551 corrCoeffs.removeAll();
00552 paramList.sort(reverse);
00553 corrCoeffs.add(paramList);
00554 }
00555 if (!In) {
00556 RooArgSet trivSet;
00557 TIterator* iter = corrCoeffs.createIterator();
00558 RooRealVar *theCorrCoef(0);
00559 while(theCorrCoef=(RooRealVar*)iter->Next()) {
00560 Double_t theVal=theCorrCoef->getVal();
00561 if ((fabs(theVal)>1)||(fabs(theVal)<1e-4))
00562 trivSet.add(*theCorrCoef);
00563 }
00564 delete iter;
00565 corrCoeffs.remove(trivSet);
00566 }
00567 params.add(corrCoeffs);
00568
00569 RooStringVar configFile("configFile", "configFile", _configFile);
00570 params.add(configFile);
00571 RooStringVar DataInputSec("DataInputSec", "DataInputSec",
00572 getDatasets()->getVarSec());
00573 params.add(DataInputSec);
00574 RooStringVar MasterPdfSec("MasterPdfSec", "MasterPdfSec", getMasterSec());
00575 params.add(MasterPdfSec);
00576 RooStringVar ActionSec("ActionSec", "ActionSec", _runSec);
00577 params.add(ActionSec);
00578
00579 if (!In) params.writeToFile(paramFile);
00580 else {
00581 Bool_t fail=params.readFromFile(paramFile);
00582 if (fail) {
00583 cout<<"Fail to read in from "<<paramFile<<endl;
00584 exit(-1);
00585 }
00586 cout<<"Info of params you just read in"<<endl
00587 <<"Config file : "<<configFile.getVal()<<endl
00588 <<"Data input Sec: "<<DataInputSec.getVal()<<endl
00589 <<"Master pdf Sec: "<<MasterPdfSec.getVal()<<endl
00590 <<"Action Sec: "<<ActionSec.getVal()<<endl;
00591 }
00592 }
00593
00599 void rarMLFitter::paramFileI(RooArgSet params, TString paramFile)
00600 {
00601 paramFileIO(params, paramFile);
00602 }
00603
00609 void rarMLFitter::paramFileO(RooArgSet params, TString paramFile)
00610 {
00611 paramFileIO(params, paramFile, kFALSE);
00612 }
00613
00618 void rarMLFitter::chkBlind(TString dsName)
00619 {
00620 if (_datasets->isBlind(dsName)) {
00621 cout<<endl<<" Dataset \""<<dsName<<"\" is still blind."<<endl
00622 <<" (any action other than pdfFit or toyStudy on the dataset"
00623 <<" requires its unblind)"<<endl
00624 <<" (unblinding dataset does not necessarily"
00625 <<" unblind your results,)"<<endl
00626 <<" (as long as you use RooUnblindPrecision or RooUnblindOffset"
00627 <<" and the state is set to blind,)"<<endl
00628 <<" (those variables they blind still remain blind.)"<<endl
00629 <<" If you really want to UNBLIND "<<dsName<<","<<endl
00630 <<" please specify in dsi section ["
00631 <<_datasets->getVarSec()<<"]:"<<endl<<endl
00632 <<" ub_"<<dsName<<" = [<OldIDsIfAny>] "<<_datasets->ubStr(dsName)
00633 <<endl<<endl<<" and re-run the job"<<endl;
00634
00635 }
00636 return;
00637 }
00638
00648 RooFitResult *rarMLFitter::doMLFit(RooDataSet *mlFitData,
00649 TString opt, ostream &o, Int_t ncpus)
00650 {
00651 cout<<endl<<" In rarMLFitter doMLFit for "<<GetName()<< " Options: " << opt
00652 << " using " << ncpus << " CPUs (if available)." << endl;
00653 RooFitResult *fitResult(0);
00654
00655 assert(_thePdf);
00656
00657
00658 Bool_t mlFitExtended = opt.Contains("e");
00659 Bool_t mlFitMinos = opt.Contains("m");
00660 Bool_t mlFitHesse = opt.Contains("h");
00661 Bool_t mlFitVerbose = !opt.Contains("q");
00662 Bool_t mlFitSave = opt.Contains("r");
00663
00664
00665
00666 fitResult=_thePdf->fitTo(*mlFitData, ConditionalObservables(_conditionalObs),
00667 Save(mlFitSave), Extended(mlFitExtended),
00668 Verbose(mlFitVerbose), Hesse(mlFitHesse),
00669 Minos(mlFitMinos), NumCPU(ncpus));
00670
00671
00672 fitResult->globalCorr();
00673
00674 saveCorrCoeffs(fitResult);
00675
00676 o<<endl<<"The mlFit results (wrt "<<mlFitData->GetName()
00677 <<" with " << mlFitData->numEntries() <<" events):"<<endl;
00678 #ifndef USENEWROOT
00679 fitResult->printToStream(o, RooPrintable::Verbose);
00680 #else
00681 Int_t contents(0);
00682 fitResult->printMultiline(o, contents, kTRUE);
00683 #endif
00684
00685 getSplitCoeffValues(getSpecialSet("cat1Set"), "Frac", o);
00686 getSplitCoeffValues(getSpecialSet("asymSet"), "Asym", o);
00687
00688 return fitResult;
00689 }
00690
00699 RooFitResult *rarMLFitter::doTheFit(RooAbsPdf *pdf, RooDataSet *fitData, TString fitOptions, Int_t ncpus)
00700 {
00701 cout<<endl<<"In rarMLFitter doTheFit for "<<GetName()<< " Options: " << fitOptions
00702 << " using " << ncpus << " CPUs (if available)." << endl;
00703
00704
00705 assert(pdf);
00706
00707
00708 Bool_t fitExtended = fitOptions.Contains("e");
00709 Bool_t fitMinos = fitOptions.Contains("m");
00710 Bool_t fitHesse = fitOptions.Contains("h");
00711 Bool_t fitVerbose = !fitOptions.Contains("q");
00712 Bool_t fitSave = fitOptions.Contains("r");
00713
00714 TStopwatch timer;
00715 timer.Start();
00716
00717 RooFitResult *fitResult=_thePdf->fitTo(*fitData, ConditionalObservables(_conditionalObs),
00718 Save(fitSave), Extended(fitExtended),
00719 Verbose(fitVerbose), Hesse(fitHesse),
00720 Minos(fitMinos), NumCPU(ncpus));
00721 timer.Stop();
00722
00723 cout<<endl<<"The doTheFit RealTime= " << timer.RealTime() << " CpuTime= "<< timer.CpuTime() << endl;
00724
00725 return fitResult;
00726 }
00727
00736 void rarMLFitter::doSignf(RooDataSet *mlFitData, TString signfStr,
00737 RooFitResult *fitResult, RooArgSet fullParams,
00738 ostream &o)
00739 {
00740 cout<<endl<<" In rarMLFitter doSignf for "<<GetName()<<endl;
00741
00742
00743 string fParamSStr0;
00744 writeToStr(fullParams, fParamSStr0);
00745
00746 TString signfFitOpt="qemhr";
00747
00748 Bool_t signfFitExtended = signfFitOpt.Contains("e");
00749 Bool_t signfFitMinos = signfFitOpt.Contains("m");
00750 Bool_t signfFitHesse = signfFitOpt.Contains("h");
00751 Bool_t signfFitVerbose = !signfFitOpt.Contains("q");
00752 Bool_t signfFitSave = signfFitOpt.Contains("r");
00753
00754
00755 if (!fitResult) {
00756 cout<<"No fit results from mlFit"<<endl;
00757 exit(-1);
00758 return;
00759 }
00760 Double_t nll=2*fitResult->minNll();
00761 o<<endl;
00762 rarStrParser signfStrParser=signfStr;
00763 while (signfStrParser.nArgs()>0) {
00764 TString paramName=signfStrParser[0];
00765 signfStrParser.Remove();
00766 RooRealVar *theParam=(RooRealVar*)fullParams.find(paramName);
00767 if (!theParam) continue;
00768
00769 readFromStr(fullParams, fParamSStr0);
00770 Double_t zSignfVal(0);
00771 if ((signfStrParser.nArgs()>0)&&(isNumber(signfStrParser[0]))) {
00772 zSignfVal = atof(signfStrParser[0]);
00773 signfStrParser.Remove();
00774 }
00775 Double_t sVal=theParam->getVal();
00776 theParam->setVal(zSignfVal);
00777 theParam->setConstant();
00778
00779
00780
00781
00782 RooFitResult *theResult=
00783 _thePdf->fitTo(*mlFitData, ConditionalObservables(_conditionalObs),
00784 Save(signfFitSave), Extended(signfFitExtended),
00785 Verbose(signfFitVerbose), Hesse(signfFitHesse),
00786 Minos(signfFitMinos));
00787
00788
00789 o<<" Signf. of "<<theParam->GetName()<<" being "<<sVal
00790 <<" wrt "<<zSignfVal<<" is "<<sqrt(2*theResult->minNll()-nll)
00791 <<" (sigma)"<<endl;
00792 }
00793
00794
00795 readFromStr(fullParams, fParamSStr0);
00796 return;
00797 }
00798
00807 void rarMLFitter::doSysStudy(RooDataSet *mlFitData, TString paramsStr,
00808 TString varsStr, RooArgSet fullParams,
00809 ostream &o)
00810 {
00811 cout<<endl<<" In rarMLFitter doSysStudy for "<<GetName()<<endl;
00812
00813
00814 string fParamSStr0;
00815 writeToStr(fullParams, fParamSStr0);
00816
00817 Int_t paramNameLen=12;
00818 {
00819 rarStrParser paramsStrParser=paramsStr;
00820 for (Int_t i=0; i<paramsStrParser.nArgs(); i++) {
00821 if (paramNameLen<paramsStrParser[i].Length())
00822 paramNameLen=paramsStrParser[i].Length();
00823 }
00824 }
00825
00826 RooArgSet studyVars;
00827 rarStrParser varsStrParser=varsStr;
00828 while (varsStrParser.nArgs()>0) {
00829 TString varName=varsStrParser[0];
00830 varsStrParser.Remove();
00831 RooRealVar *theVar=(RooRealVar*)fullParams.find(varName);
00832 if (theVar) studyVars.add(*theVar);
00833 }
00834 if (studyVars.getSize()<1) {
00835 cout<<" No vars in postMLSysVars found!"<<endl;
00836 return;
00837 }
00838 cout<<" postMLSysVars:"<<endl;
00839 studyVars.Print("v");
00840
00841 TString sysFitOpt="qemhr";
00842
00843
00844
00845
00846 Bool_t sysFitExtended = sysFitOpt.Contains("e");
00847 Bool_t sysFitMinos = sysFitOpt.Contains("m");
00848 Bool_t sysFitHesse = sysFitOpt.Contains("h");
00849 Bool_t sysFitVerbose = !sysFitOpt.Contains("q");
00850 Bool_t sysFitSave = sysFitOpt.Contains("r");
00851
00852 _thePdf->fitTo(*mlFitData, ConditionalObservables(_conditionalObs),
00853 Save(sysFitSave), Extended(sysFitExtended),
00854 Verbose(sysFitVerbose), Hesse(sysFitHesse),
00855 Minos(sysFitMinos));
00856
00857
00858 string fParamSStr;
00859 writeToStr(fullParams, fParamSStr);
00860
00861 RooArgSet *cStudyVars=(RooArgSet*)studyVars.snapshot(kTRUE);
00862
00863 rarStrParser paramsStrParser=paramsStr;
00864
00865 Double_t defVU(1);
00866 if (paramsStrParser.nArgs()>0) {
00867 if (isNumber(paramsStrParser[0])) {
00868 defVU=fabs(atof(paramsStrParser[0]));
00869 paramsStrParser.Remove();
00870 }
00871 }
00872 Int_t nSysParams(0);
00873
00874 TString vParams="";
00875 TArrayD pV(1), mV(1);
00876 TArrayD pArray(studyVars.getSize()), mArray(studyVars.getSize());
00877 TArrayD aArray(studyVars.getSize());
00878
00879 while (paramsStrParser.nArgs()>0) {
00880 TString theParamName=paramsStrParser[0];
00881 paramsStrParser.Remove();
00882 RooArgSet theParamSet;
00883 TIterator* iter = fullParams.createIterator();
00884 RooRealVar *theParam(0);
00885 Bool_t foundvParam=kFALSE;
00886 while(theParam=(RooRealVar*)iter->Next()) {
00887 TString theName=theParam->GetName();
00888 if ((theName==theParamName)||(theName.BeginsWith(theParamName+"_"))) {
00889 theParamSet.add(*theParam);
00890 if (!foundvParam) {
00891 foundvParam=kTRUE;
00892 vParams=vParams+" \""+theParamName+"\"";
00893 }
00894 }
00895 }
00896 delete iter;
00897 if (theParamSet.getSize()<1) continue;
00898 nSysParams++;
00899
00900 pV.Set(nSysParams);
00901 mV.Set(nSysParams);
00902 pArray.Set(nSysParams*studyVars.getSize());
00903 mArray.Set(nSysParams*studyVars.getSize());
00904 aArray.Set(nSysParams*studyVars.getSize());
00905
00906 cout<<" SysStudy for "<<theParamName<<endl;
00907 theParamSet.Print("v");
00908
00909 Bool_t useErr=kTRUE;
00910 pV[nSysParams-1]=defVU;
00911
00912 if (paramsStrParser.nArgs()>0) {
00913 if (isNumber(paramsStrParser[0])) {
00914 TString myVStr=paramsStrParser[0];
00915 paramsStrParser.Remove();
00916 pV[nSysParams-1]=fabs(atof(myVStr));
00917 if (myVStr.EndsWith("V")||(myVStr.EndsWith("v"))) useErr=kFALSE;
00918 }
00919 }
00920
00921 readFromStr(fullParams, fParamSStr);
00922
00923 setVariation(theParamSet, pV[nSysParams-1], useErr, kTRUE);
00924
00925
00926 _thePdf->fitTo(*mlFitData, ConditionalObservables(_conditionalObs),
00927 Save(sysFitSave), Extended(sysFitExtended),
00928 Verbose(sysFitVerbose), Hesse(sysFitHesse),
00929 Minos(sysFitMinos));
00930
00931
00932
00933 calSysErrors(nSysParams-1, *cStudyVars, studyVars, pArray);
00934
00935
00936 mV[nSysParams-1]=pV[nSysParams-1];
00937 if (paramsStrParser.nArgs()>0) {
00938 if (isNumber(paramsStrParser[0])) {
00939 TString myVStr=paramsStrParser[0];
00940 paramsStrParser.Remove();
00941 mV[nSysParams-1]=fabs(atof(myVStr));
00942 if (myVStr.EndsWith("V")||(myVStr.EndsWith("v"))) useErr=kFALSE;
00943 }
00944 }
00945
00946 readFromStr(fullParams, fParamSStr);
00947
00948 setVariation(theParamSet, mV[nSysParams-1], useErr, kFALSE);
00949
00950 _thePdf->fitTo(*mlFitData, ConditionalObservables(_conditionalObs),
00951 Save(sysFitSave), Extended(sysFitExtended),
00952 Verbose(sysFitVerbose), Hesse(sysFitHesse),
00953 Minos(sysFitMinos));
00954
00955
00956 calSysErrors(nSysParams-1, *cStudyVars, studyVars, mArray);
00957 }
00958
00959 TMatrixD corrM1=getCorrMatrix(nSysParams, vParams, kTRUE);
00960 TMatrixD corrM=getCorrMatrix(nSysParams, vParams);
00961
00962 o<<endl<<" Systematic Error Table:"<<endl;
00963 o<<setw(paramNameLen+11)<<" ";
00964
00965 RooArgList studyVarList(studyVars);
00966 for (Int_t i=0; i<studyVarList.getSize(); i++) {
00967 o<<setw(40)<<studyVarList[i].GetName();
00968 }
00969
00970 o<<setw(15)<<"corr matrix:";
00971 rarStrParser vParamsParser=vParams;
00972 for (Int_t i=0; i<nSysParams; i++) {
00973 o<<setw(paramNameLen+1)<<vParamsParser[i];
00974 }
00975 o<<endl;
00976
00977 avgSysErrors(o, vParams, paramNameLen, pV,pArray, mV,mArray, aArray, corrM);
00978
00979 outSysErrors(o, "(w/o corr):", paramNameLen, corrM1, aArray);
00980 outSysErrors(o, "(w/ corr):", paramNameLen, corrM, aArray);
00981
00982 readFromStr(fullParams, fParamSStr0);
00983 return;
00984 }
00985
00991 void rarMLFitter::setVariation(RooArgSet theParams, Double_t myV,
00992 Bool_t useErr, Bool_t isPlus)
00993 {
00994 TIterator* iter = theParams.createIterator();
00995 RooRealVar *theParam(0);
00996 while(theParam=(RooRealVar*)iter->Next()) {
00997 Double_t myVariant=myV;
00998 if (useErr) {
00999 Double_t myErr(0);
01000 if (theParam->hasAsymError()) {
01001 if (isPlus) myErr=theParam->getAsymErrorHi();
01002 else myErr=theParam->getAsymErrorLo();
01003 }
01004 if ((0==myErr)&&(theParam->hasError()))
01005 myErr=theParam->getError();
01006 if (0==myErr) {
01007 cout<<" No err specified for "<<theParam->GetName()<<endl;
01008 continue;
01009 }
01010 myVariant=myV*myErr;
01011 }
01012 myVariant=fabs(myVariant);
01013 if (!isPlus) myVariant=-myVariant;
01014 cout<<" Variation for "<<theParam->GetName()<<": "<<myVariant<<endl;
01015 theParam->setVal(theParam->getVal()+myVariant);
01016 }
01017 delete iter;
01018 return;
01019 }
01020
01026 void rarMLFitter::calSysErrors(Int_t iParam, RooArgSet &cstudyVars,
01027 RooArgSet &studyVars, TArrayD &eArray)
01028 {
01029 TIterator* iter = studyVars.createIterator();
01030 RooRealVar *theVar(0), *theIVar(0);
01031 Int_t vIdx(iParam*studyVars.getSize());
01032 while(theVar=(RooRealVar*)iter->Next()) {
01033 theIVar=(RooRealVar*)cstudyVars.find(theVar->GetName());
01034 if (!theIVar) {
01035 cout<<" Can not find "<<theVar->GetName()<<endl;
01036 exit(-1);
01037 }
01038 eArray[vIdx++]=theVar->getVal()-theIVar->getVal();
01039 }
01040
01041 delete iter;
01042 return;
01043 }
01044
01055 void rarMLFitter::avgSysErrors(ostream &o, TString vParams, Int_t pnLen,
01056 TArrayD &pV, TArrayD &pArray,
01057 TArrayD &mV, TArrayD &mArray,
01058 TArrayD &aArray, TMatrixD corrM)
01059 {
01060 Int_t nParam=pV.GetSize();
01061 rarStrParser vParamsParser=vParams;
01062 if (nParam!=vParamsParser.nArgs()) {
01063 cout<<vParams<<endl
01064 <<"does not match # of params: "<<nParam<<endl;
01065 o<<vParams<<endl
01066 <<"does not match # of params: "<<nParam<<endl;
01067 }
01068 Int_t nVar=pArray.GetSize()/nParam;
01069 for (Int_t j=0; j<nParam; j++) {
01070 o.unsetf(ios_base::scientific);
01071 o.setf(ios_base::showpos);
01072 o.precision(4);
01073 o<<setw(pnLen)<<vParamsParser[j]<<setw(5)<<pV[j]<<setw(5)<<-mV[j]<<":";
01074 o.setf(ios_base::scientific);
01075 for (Int_t i=0; i<nVar; i++) {
01076
01077 o<<setw(16)<<pArray[j*nVar+i]
01078 <<setw(12)<<mArray[j*nVar+i];
01079
01080 Double_t aError(0);
01081 if (pV[j]>0) aError+=pArray[j*nVar+i];
01082 if (mV[j]>0) aError-=mArray[j*nVar+i];
01083 if ((pV[j]>0)&&(mV[j]>0)) {
01084 aError/=2.;
01085 if (pArray[j*nVar+i]*mArray[j*nVar+i]>0) {
01086 cout<<" W A R N I N G !"<<endl
01087 <<" + and - variation have the same sign!"<<endl
01088 <<" Using the larger one only!"<<endl;
01089 aError=pArray[j*nVar+i];
01090 if (fabs(aError)<fabs(mArray[j*nVar+i]))
01091 aError=mArray[j*nVar+i];
01092 }
01093 }
01094 aArray[j*nVar+i]=aError;
01095 o.unsetf(ios_base::showpos);
01096 o<<setw(12)<<aError;
01097 }
01098
01099 o<<setw(15)<<" ";
01100 for (Int_t i=0; i<nParam; i++) {
01101 o<<setw(pnLen+1)<<corrM(j, i);
01102 }
01103 o<<endl;
01104 }
01105
01106 return;
01107 }
01108
01115 void rarMLFitter::outSysErrors(ostream &o, TString rowName, Int_t pnLen,
01116 TMatrixD corrM, TArrayD &aArray)
01117 {
01118
01119 o<<setw(pnLen)<<rowName<<setw(11)<<" ";
01120 o.setf(ios_base::scientific);
01121
01122 Int_t nSysParams=corrM.GetNcols();
01123 Int_t nSysVars=aArray.GetSize()/nSysParams;
01124 for (Int_t i=0; i<nSysVars; i++) {
01125
01126 TMatrixD errM(nSysParams,1);
01127 for (Int_t j=0; j<nSysParams; j++) {
01128 errM(j,0)=aArray[j*nSysVars+i];
01129 }
01130
01131 TMatrixD errMT(1,nSysParams);
01132 errMT.Transpose(errM);
01133
01134 TMatrixD errorSq(errMT*corrM*errM);
01135 o<<setw(40)<<sqrt(errorSq(0,0));
01136
01137 }
01138 o<<endl;
01139
01140 }
01141
01146 TMatrixD rarMLFitter::getCorrMatrix(Int_t nSysParams, TString vParams,
01147 Bool_t diagOnly)
01148 {
01149 TMatrixD retM(nSysParams, nSysParams);
01150 retM.Zero();
01151 rarStrParser vParamsParser=vParams;
01152 if (vParamsParser.nArgs()!=nSysParams) {
01153 cout<<" Names and number of names ("<<nSysParams<<") do not match"<<endl
01154 <<vParams<<endl;
01155 return retM;
01156 }
01157 for (Int_t i=0; i<nSysParams; i++) {
01158 for (Int_t j=0; j<nSysParams; j++) {
01159 if (i==j) {
01160 retM(i,j)=1;
01161 continue;
01162 }
01163 if (diagOnly) continue;
01164 retM(i,j)=getCorrCoeff(vParamsParser[i], vParamsParser[j]);
01165 }
01166 }
01167 return retM;
01168 }
01169
01177 Double_t rarMLFitter::doGOFChisq(RooDataSet *mlFitData, ostream &o,
01178 TList *plotList)
01179 {
01180
01181 RooArgSet GOFObsSet;
01182 addProtVars("postMLGOFObs", GOFObsSet);
01183 GOFObsSet.add(*_thePdf->getDependents(mlFitData), kTRUE);
01184 if (GOFObsSet.getSize()<=0) return 0;
01185
01186 string obsStr;
01187 writeToStr(GOFObsSet, obsStr);
01188
01189 RooArgList GOFObsList(GOFObsSet);
01190 for (Int_t i=0; i<GOFObsList.getSize(); i++) {
01191 RooRealVar *theVar=dynamic_cast<RooRealVar*>(GOFObsList.at(i));
01192 if (!theVar) continue;
01193 Double_t min, max;
01194 getRange(theVar, "plotRange_", min, max, _runSec);
01195 Int_t nBins=atoi
01196 (readConfStr(Form("plotBins_%s",theVar->GetName()),"-1",_runSec));
01197 if (nBins>0) theVar->setBins(nBins);
01198 }
01199 if (plotList) {
01200 cout<<"GOFObsSet:"<<endl;
01201 GOFObsSet.Print("v");
01202 }
01203
01204
01205
01206 RooDataHist dHist("dHist", "dHist", GOFObsSet, *mlFitData);
01207 if (plotList) {
01208
01209 dHist.dump2();
01210 }
01211
01212
01213
01214
01215
01216 RooSuperCategory *theSuperCat(0);
01217 RooArgSet inputCats;
01218 if (getControlBit("SimFit")) {
01219 theSuperCat=(RooSuperCategory*) &((RooSimultaneous*)_thePdf)->indexCat();
01220 inputCats.add(theSuperCat->inputCatList());
01221 }
01222
01223 RooDataHist pHist("pHist", "pHist", GOFObsSet);
01224 Int_t nBins=pHist.numEntries();
01225 Double_t sum(0);
01226 for (Int_t i=0; i<nBins; i++) {
01227 const RooArgSet *theBinSet=pHist.get(i);
01228
01229 RooArgList theBinList(*theBinSet);
01230 for (Int_t j=0; j<theBinList.getSize(); j++) {
01231 RooRealVar *theBin=dynamic_cast<RooRealVar*>(theBinList.at(j));
01232 if (!theBin) {
01233 RooAbsCategory *theCat=dynamic_cast<RooAbsCategory*>
01234 (theBinList.at(j));
01235 if (!theCat) continue;
01236 RooAbsCategoryLValue *theInCat=(RooAbsCategoryLValue*)
01237 inputCats.find(theCat->GetName());
01238 if (!theInCat) continue;
01239 theInCat->setLabel(theCat->getLabel());
01240 continue;
01241 }
01242 RooRealVar *theVar=dynamic_cast<RooRealVar*>(GOFObsList.at(j));
01243 assert(theVar);
01244 Double_t hbSize=(theBin->getMax()-theBin->getMin())/theBin->getBins()/2.;
01245 theVar->setRange("subrange",
01246 theBin->getVal()-hbSize, theBin->getVal()+hbSize);
01247 theVar->setVal(theBin->getVal());
01248
01249
01250 }
01251
01252 RooAbsPdf *theIntPdf(_thePdf);
01253 if (theSuperCat) {
01254 theIntPdf=((RooSimultaneous*)_thePdf)->getPdf(theSuperCat->getLabel());
01255
01256
01257 }
01258
01259 RooAbsReal *pdfIntegral(0);
01260 if (theIntPdf) pdfIntegral=
01261 theIntPdf->createIntegral(GOFObsSet, GOFObsSet, "subrange");
01262
01263 Double_t thisIntegral(0);
01264 RooArgSet nullDS;
01265 if (pdfIntegral) thisIntegral=
01266 theIntPdf->expectedEvents(&nullDS)*pdfIntegral->getVal();
01267 if (theSuperCat) thisIntegral/=theSuperCat->numTypes();
01268 sum+=thisIntegral;
01269
01270
01271 pHist.set(*theBinSet, thisIntegral);
01272
01273 }
01274 if (plotList) {
01275 cout<<"sum="<<sum<<endl;
01276 pHist.dump2();
01277 }
01278
01279 TH1F *dataHist(0);
01280 TH1F *pdfHist(0);
01281 if (plotList) {
01282 dataHist=new TH1F("dataHist_GOF", "dataHist_GOF", nBins, 0, 1);
01283 pdfHist=new TH1F("pdfHist_GOF", "pdfHist_GOF", nBins, 0, 1);
01284 plotList->Add(dataHist);
01285 plotList->Add(pdfHist);
01286 }
01287
01288 Double_t chisq(0);
01289
01290 for (Int_t i=0; i<nBins; i++) {
01291 const RooArgSet *theBinSet=dHist.get(i);
01292
01293 Double_t dw=dHist.weight(*theBinSet,0);
01294 Double_t ym,yp;
01295 RooHistError::instance().getPoissonInterval(Int_t(dw+0.5),ym,yp,1);
01296 Double_t del=dw-ym;
01297 Double_t deh=yp-dw;
01298 Double_t pw=pHist.weight(*theBinSet,0);
01299 if (plotList) {
01300 dataHist->SetBinContent(i, dw);
01301 pdfHist->SetBinContent(i, pw);
01302 }
01303 Double_t pwmdw=pw-dw;
01304 Double_t err=(pwmdw>0)?deh:del;
01305 if (0==err) continue;
01306 chisq+=pwmdw*pwmdw/(err*err);
01307
01308 if (plotList)
01309 cout<<"err="<<err<<" dw="<<dw<<" pw="<<pw<<endl;
01310
01311 }
01312 o<<"GOF chisq: "<<chisq<<endl;
01313
01314
01315 readFromStr(GOFObsSet, obsStr);
01316 return chisq;
01317 }
01318
01334 RooDataSet *rarMLFitter::doToyStudy(RooArgSet fullParams)
01335 {
01336 rarToyList toylist;
01337
01338 cout<<endl<<" In rarMLFitter doToyStudy for "<<GetName()<<endl;
01339
01340 TString temp=readConfStr("toyVerbose", "yes", _runSec);
01341 Bool_t toyVerbose = kFALSE;
01342 if (temp == "yes") {toyVerbose = kTRUE;}
01343
01344 RooDataSet *toyResults(0);
01345
01346 TString protDataStr=readConfStr("protToyData", "", _runSec);
01347 if (""==protDataStr) protDataStr=readConfStr("protDatasets", "", _runSec);
01348 rarStrParser protDataStrParser=protDataStr;
01349
01350 if (protDataStrParser.nArgs()>0) {
01351 TString mpdsName=protDataStrParser[0];
01352 protDataStrParser.Remove();
01353 _protDataset=_datasets->getData(mpdsName);
01354 if (!_protDataset) {
01355 cout<<" Toy Error: Can not find dataset named "<<mpdsName<<endl
01356 <<" for your prototype dataset."<<endl
01357 <<" Please make sure you have the right Dataset name"<<endl;
01358 exit(-1);
01359 }
01360 } else if (_theData) _protDataset=_theData;
01361 else {
01362 cout<<" Toy Error: No prototype datasets defined!!"<<endl
01363 <<" You could have problem with your toy study"<<endl
01364 <<" See config `protDatasets' in online doc for more info"<<endl;
01365 exit(-1);
01366 }
01367 cout<<" Toy Using "<<_protDataset->GetName()<<" as master protDataset"<<endl;
01368
01369
01370 TString preToyRandParams="";
01371
01372 RooArgSet paramRandSet;
01373 if (_theToyParamGen) {
01374 RooArgSet paramNobs(*_thePdf->getParameters(_protDataset));
01375 paramNobs.add(*_thePdf->getDependents(_protDataset));
01376 RooArgList allOtherParams(*_theToyParamGen->getParameters(paramNobs));
01377 for (Int_t i=0; i<allOtherParams.getSize(); i++) {
01378 RooRealVar *theParam=(RooRealVar*)allOtherParams.at(i);
01379 if (theParam->hasMin()&&theParam->hasMax()) paramRandSet.add(*theParam);
01380 }
01381
01382 preToyRandParams=readConfStr("preToyRandParams", "", _runSec);
01383 }
01384 rarStrParser preToyRandParser=preToyRandParams;
01385 if ((preToyRandParser.nArgs()>0)&&(preToyRandParser[0]=="no"))
01386 preToyRandParser="";
01387
01388
01389 Int_t toyNexp=atoi(readConfStr("toyNexp", "1", _runSec));
01390 if (_toyNexp>0) toyNexp=_toyNexp;
01391
01392 rarStrParser toyNevtParser=readConfStr("toyNevt", "0 fixed", _runSec);
01393 Double_t toyNevt=atoi(toyNevtParser[0]);
01394 #ifndef USENEWROOT
01395 if (toyNevt<=0) toyNevt=_protDataset->numEntries(kTRUE);
01396 #else
01397 if (toyNevt<=0) toyNevt=_protDataset->numEntries();
01398 #endif
01399 toylist.setDataset(_protDataset->GetName(), toyNevt);
01400
01401 Bool_t extendedGen=
01402 ("floated"==toyNevtParser[1])||
01403 ("extended"==toyNevtParser[1])||
01404 ("notfixed"==toyNevtParser[1]);
01405 cout<<"Toy Total Number of events: "<<toyNevt << " to be generated ";
01406 if (extendedGen) cout<<"with";
01407 else cout<<"without";
01408 cout<<" Poisson fluctuation"<<endl;
01409
01410 RooArgSet compCoeffSet("Comp Coeff Set");
01411 RooArgSet coeffParamSet("Comp Coeff Parameter Set");
01412 Int_t nComp=1;
01413 if (getControlBit("SimFit")) nComp=_nComp;
01414 for (Int_t i=0; i<nComp; i++) {
01415 rarAdd *compPdf=(rarAdd*)_pdfList.At(i);
01416 RooArgList thisCompCoeffList(compPdf->getCoeffList());
01417
01418 compCoeffSet.add(thisCompCoeffList);
01419 cout << "Toy Size thisCompCoeffList " << thisCompCoeffList.getSize() << endl;
01420 thisCompCoeffList.Print("v");
01421 for (Int_t j=0; j<thisCompCoeffList.getSize(); j++) {
01422
01423
01424 RooArgList thisCompCoeffParams(*(thisCompCoeffList[j].getParameters(*_fullObs)));
01425
01426
01427 if (thisCompCoeffParams.getSize() == 0) {
01428 coeffParamSet.add(thisCompCoeffList[j]);
01429 } else {
01430 for (Int_t k=0; k<thisCompCoeffParams.getSize(); k++) {
01431 TString thisParamName=thisCompCoeffParams[k].GetName();
01432
01433 if (isFracName(thisParamName)) continue;
01434 coeffParamSet.add(thisCompCoeffParams[k]);
01435 }
01436 }
01437 }
01438 }
01439 Int_t nCoeff=compCoeffSet.getSize();
01440 Int_t nCoeffParam=coeffParamSet.getSize();
01441 RooArgList compCoeffList(compCoeffSet);
01442 RooArgList coeffParamList(coeffParamSet);
01443
01444 if (toyVerbose) {
01445 cout << "Toy: List of coefficients: " << nCoeff << endl;
01446 compCoeffList.Print("v");
01447 cout << "Toy: List of parameters: " << nCoeffParam << endl;
01448 coeffParamList.Print("v");
01449 }
01450
01451
01452 RooArgSet evtTypeVarSet("Event Type Var Set");
01453 RooArgSet evtTypeStrSet("Event Type Str Set");
01454 for (Int_t i=0; i<nCoeffParam; i++) {
01455 RooAbsReal *theCoeff=dynamic_cast<RooAbsReal*>(coeffParamList.at(i));
01456 if (!theCoeff) {
01457 cout<<endl<<"Toy W A R N I N G ! ! !"<<endl
01458 <<" "<<coeffParamList[i].GetName()
01459 <<" is not a (sub-)class of RooAbsReal."<<endl
01460 <<" If it is a blind category var,"
01461 <<" please set it to unblind for toy study."<<endl
01462 <<" ie, by adding in section ["<<_runSec<<"]"<<endl
01463 <<" "<<coeffParamList[i].GetName()<<" = 0"<<endl
01464 <<endl;
01465 continue;
01466 }
01467
01468 TString coeffName = theCoeff->GetName();
01469 toylist.setInitial(coeffName, theCoeff->getVal());
01470
01471 rarStrParser toySrcParser=
01472 readConfStr(Form("toySrc_%s", theCoeff->GetName()),
01473 Form("pdf %f",theCoeff->getVal()), _runSec);
01474 Double_t coeffParamVal(0);
01475 Int_t nSrcType=toySrcParser.nArgs()/2;
01476 for (Int_t j=0; j<nSrcType; j++) {
01477 TString evtSrc=toySrcParser[0];
01478 toySrcParser.Remove();
01479 TString evtSrcStr=toySrcParser[0];
01480 toySrcParser.Remove();
01481 Double_t evtSrcVal=getFormulaVal(evtSrcStr);
01482 coeffParamVal+=evtSrcVal;
01483 TString evtTypeName=
01484 Form("%s \"%s\"", coeffParamList[i].GetName(), evtSrc.Data());
01485 RooRealVar *evtTypeVar=(RooRealVar*)evtTypeVarSet.find(evtTypeName);
01486 toylist.setMethod(coeffName, evtSrc.Data());
01487
01488 if (!evtTypeVar) {
01489 evtTypeVar=new RooRealVar(evtTypeName, evtTypeName, 0);
01490 evtTypeVarSet.addOwned(*evtTypeVar);
01491 }
01492 evtTypeVar->setVal(evtTypeVar->getVal()+evtSrcVal);
01493
01494 RooStringVar *evtTypeStr=(RooStringVar*)evtTypeStrSet.find(evtTypeName);
01495 if (!evtTypeStr) {
01496 evtTypeStr=new RooStringVar(evtTypeName, evtTypeName, "");
01497 evtTypeStrSet.addOwned(*evtTypeStr);
01498 }
01499 evtTypeStr->
01500 setVal(Form("%s \"%s\"",evtTypeStr->getVal(),evtSrcStr.Data()));
01501 }
01502 RooRealVar *coeffParamVar=(RooRealVar*)coeffParamList.at(i);
01503
01504 if (!coeffParamVar->inRange(coeffParamVal, "range")) {
01505 cout<<coeffParamVal<<" not within the range of "
01506 <<coeffParamVar->GetName()<<": ("
01507 <<coeffParamVar->getMin()<<","
01508 <<coeffParamVar->getMax()<<")"<<endl;
01509 exit(-1);
01510 }
01511 coeffParamVar->setVal(coeffParamVal);
01512
01513 toylist.setRequested(coeffName, coeffParamVal);
01514 }
01515
01516 if (toyVerbose) {
01517 cout << "Toy: List of generating methods: " << nCoeff << endl;
01518 evtTypeVarSet.Print("v");
01519 evtTypeStrSet.Print("v");
01520 }
01521
01522 RooArgList evtTypeVarList(evtTypeVarSet);
01523 RooArgList evtTypeStrList(evtTypeStrSet);
01524 Int_t nEvtType=evtTypeVarList.getSize();
01525 {
01526
01527
01528 TString toyAdjustMethod=readConfStr("toyAdjustMethod", "Largest", _runSec);
01529 toyAdjustMethod.ToLower();
01530 if (toyAdjustMethod.Contains("largest") ||
01531 toyAdjustMethod.Contains("total") ||
01532 toyAdjustMethod.Contains("smallest") ||
01533 toyAdjustMethod.Contains("prorata")) {
01534
01535
01536 } else {
01537 cout << "Toy toyAdjustMethod: expected token [Largest|Total|Smallest|Prorata]"
01538 << " found : " << toyAdjustMethod << "." << endl;
01539
01540 }
01541
01542
01543
01544
01545 TString zEvtName("notSet"), lEvtName("notSet");
01546 Int_t zEvtIdx(-1),lEvtIdx(-1);
01547 Double_t lEvt(0);
01548 for (Int_t i=0; i<nEvtType; i++) {
01549 RooRealVar *evtTypeVar=(RooRealVar *)evtTypeVarList.at(i);
01550 TString evtTypeVarName=evtTypeVar->GetName();
01551 Double_t evtTypeVarVal=evtTypeVar->getVal();
01552 if(0==evtTypeVarVal) {
01553 if (zEvtIdx!=-1) {
01554 cout<<"Toy Can not have more than one coeff=0"<<endl;
01555 exit(-1);
01556 } else {
01557
01558 zEvtName=evtTypeVarName;
01559 }
01560 }
01561 if(evtTypeVarVal>lEvt) {
01562 lEvt=evtTypeVarVal;
01563 lEvtIdx=i;
01564 lEvtName=evtTypeVarName;
01565 }
01566 }
01567 if (toyVerbose) {
01568 if (zEvtName != "notSet") {
01569 cout << "Toy Dataset with #evts = 0 specified: " << zEvtName << endl;
01570 }
01571 cout << "Toy Dataset with largest #evts specified: " << lEvtName << endl;
01572 }
01573
01574 Bool_t zParamAdj(kFALSE), lParamAdj(kFALSE);
01575 while (1) {
01576 if (-1==zEvtIdx) break;
01577 if (!compCoeffSet.find(rarStrParser(zEvtName)[0])) {
01578 if (!compCoeffSet.find("theSim_"+rarStrParser(zEvtName)[0])) {
01579
01580 cout<<"Toy " <<zEvtName<<" is set to zero but is not adjustable"<<endl;
01581 exit(-1);
01582 break;
01583 }
01584 }
01585 zParamAdj=kTRUE;
01586 break;
01587 }
01588 while (1) {
01589 if (-1==lEvtIdx) break;
01590 if (!compCoeffSet.find(rarStrParser(lEvtName)[0])) {
01591 if (!compCoeffSet.find("theSim_"+rarStrParser(lEvtName)[0])) {
01592
01593 cout<<"Toy "<<lEvtName<<" is the largest but might not be adjustable"<<endl;
01594 }
01595 }
01596 lParamAdj=kTRUE;
01597 break;
01598 }
01599
01600 Double_t totEvt(0);
01601 for (Int_t i=0; i<nCoeff; i++) { totEvt+=((RooAbsReal&)compCoeffList[i]).getVal(); }
01602
01603 if (totEvt!=toyNevt) {
01604 Double_t nAdjust=toyNevt-totEvt;
01605 Int_t adjIdx(-1);
01606
01607 if (zParamAdj) adjIdx=zEvtIdx;
01608 else if (lParamAdj) adjIdx=lEvtIdx;
01609 else {
01610 cout << "Toy #totEvts = " << totEvt << " #toyNevt = " << toyNevt << endl;
01611 cout<<"Need to adjust number of evenst but none of the toySrc_* can be adjusted"
01612 <<endl;
01613 exit(-1);
01614 }
01615
01616
01617 RooRealVar *evtTypeVar=(RooRealVar *)evtTypeVarList.at(adjIdx);
01618 evtTypeVar->setVal(evtTypeVar->getVal()+nAdjust);
01619
01620 RooRealVar *coeffAdjust=(RooRealVar*)
01621 compCoeffSet.find(rarStrParser(evtTypeVar->GetName())[0]);
01622 if (!coeffAdjust) coeffAdjust=(RooRealVar*)
01623 coeffParamSet.find(rarStrParser(evtTypeVar->GetName())[0]);
01624 if (!coeffAdjust) {
01625 cout<<"Toy Can not find var "<<rarStrParser(evtTypeVar->GetName())[0]
01626 <<endl;
01627 exit(-1);
01628 }
01629 nAdjust+=coeffAdjust->getVal();
01630 cout<<"Toy "<<coeffAdjust->GetName()<<" needs to be adjusted from "
01631 <<coeffAdjust->getVal()<<" to "<<nAdjust<<endl;
01632 if (!coeffAdjust->inRange(nAdjust, "range")) {
01633 cout<<"Toy "<<nAdjust<<" not within the range of "<<coeffAdjust->GetName()
01634 <<": ("<<coeffAdjust->getMin()<<","
01635 <<coeffAdjust->getMax()<<")"<<endl;
01636 exit(-1);
01637 }
01638 coeffAdjust->setVal(nAdjust);
01639 {
01640 Double_t totEvt(0);
01641 for (Int_t i=0; i<nCoeff; i++) {totEvt+=((RooAbsReal&)compCoeffList[i]).getVal();}
01642
01643 if (fabs(totEvt-toyNevt)>.01) {
01644 cout<<endl<<"Toy Adjusting "<<evtTypeVar->GetName()
01645 <<" will not set totEvtGen "<<totEvt<<" to "<<toyNevt<<endl
01646 <<" If you are using blind var for yields,"<<endl
01647 <<" please make sure you set it to unblind"<<endl
01648 <<" (or temporarily do not use blind var for toy study)"<<endl
01649 <<" (search for `W A R N I N G ! ! !\' in your log file)"<<endl;
01650 exit(-1);
01651 }
01652 }
01653 }
01654
01655
01656 for (Int_t i=0; i<nCoeff; i++) {
01657 Double_t nevts = ((RooAbsReal&)compCoeffList[i]).getVal();
01658 TString coeffName = ((RooAbsReal&)compCoeffList[i]).getTitle();
01659
01660 toylist.setUsed(coeffName, nevts);
01661 }
01662
01663
01664 }
01665 toylist.print();
01666
01667
01668 TString toyFilePrefix=readConfStrCnA("toyDataFilePrefix", "default");
01669 if (("default"==toyFilePrefix)||
01670 ("no"==toyFilePrefix))
01671 toyFilePrefix="no";
01672 else {
01673 if ("yes"==toyFilePrefix) toyFilePrefix=Form("N \"%s\"", _runSec.Data());
01674 toyFilePrefix=getParamFileName("toySample", toyFilePrefix, _toyDir);
01675 toyFilePrefix+=Form(".%03d", _toyID);
01676 toyFilePrefix+=".%03d";
01677 }
01678
01679 TString genOpt="r";
01680 if (extendedGen) genOpt+="e";
01681 TString fitOpt=readConfStr("toyFitOption", "mhq", _runSec);
01682
01683 fitOpt+="er";
01684
01685
01686 Bool_t toyFitMinos = fitOpt.Contains("m");
01687 Bool_t toyFitHesse = fitOpt.Contains("h");
01688 Bool_t toyFitVerbose = !fitOpt.Contains("q");
01689
01690
01691
01692 Int_t toyFitNumCPU=atoi(readConfStr("useNumCPU", "1", getMasterSec()));
01693
01694
01695 RooArgSet toyFitMinosAS;
01696 TString toyFitMinosStr=readConfStr("toyFitMinos", "notSet", _runSec);
01697 if ("notSet"!=toyFitMinosStr) {
01698 toyFitMinosAS.add(getArgSet(toyFitMinosStr,kFALSE,&fullParams));
01699 }
01700
01701
01702 if (toyVerbose) {
01703 cout<<"Toy full obs for toy study"<<endl;
01704 _fObsSet.Print("v");
01705 }
01706
01707
01708 if (!getProtGen()) {
01709 cout<<"Toy Can not create protGen!"<<endl;
01710 exit(-1);
01711 }
01712
01713
01714 getCompCatDS(&_protDatasetsM, _protDataset);
01715
01716 const char *dummy("");
01717 for (Int_t i=0; i<_nComp; i++) {
01718 rarMLPdf *model=(rarMLPdf*)_pdfList.At(i);
01719 TList *modelPdfList=model->getPdfList();
01720 Int_t nModelComp=modelPdfList->GetSize();
01721 for (Int_t j=0; j<nModelComp; j++) {
01722 rarBasePdf *comp=(rarBasePdf*)modelPdfList->At(j);
01723 RooDataSet *theData(0);
01724 if (protDataStrParser.nArgs()>0) {
01725 theData=_datasets->getData(protDataStrParser[0]);
01726 protDataStrParser.Remove();
01727 } else theData=comp->getData(dummy);
01728 _compCat.setIndex(j);
01729 getCompCatDS(&_protDatasets, theData, &_compCat);
01730 }
01731 }
01732
01733
01734 _protGenLevel=atoi(readConfStr("protDataGenLevel", "1", _runSec));
01735 if (_protGenLevel>3) _protGenLevel=3;
01736 if (_protGenLevel<=0) _protGenLevel=0;
01737
01738 RooArgSet fullProtVars(_protDataVars);
01739 {
01740 RooArgSet protGenVars(*_theProtGen->getDependents(_protDataset));
01741 fullProtVars.add(protGenVars);
01742 if (fullProtVars.getSize()<=0) _protGenLevel=0;
01743 else {
01744 if (0==_protGenLevel) _protGenLevel=2;
01745 }
01746
01747 if (toyVerbose) {
01748 if (fullProtVars.getSize() > 0) {
01749 cout<<"Toy Full protDataVars defined:"<<endl;
01750 fullProtVars.Print("v");
01751 } else {
01752 cout<<"Toy No Full protDataVars defined."<<endl;
01753 }
01754 cout<<"Toy protGenLevel: "<<_protGenLevel<<endl;
01755 }
01756 }
01757
01758 if (!getGenerator()) {
01759 cout<<"Toy Can not create toy Gen!"<<endl;
01760 exit(-1);
01761 }
01762
01763 RooArgSet pdfNparamFSet(*_theGen, *_theProtGen);
01764 pdfNparamFSet.add(_subPdfs);
01765 RooArgSet *fCloneSet=(RooArgSet*)pdfNparamFSet.snapshot(kTRUE);
01766 RooAbsPdf *theGen=(RooAbsPdf*)fCloneSet->find(_theGen->GetName());
01767 RooAbsPdf *theProtGen=(RooAbsPdf*)fCloneSet->find(_theProtGen->GetName());
01768
01769
01770 RooDataSet *protData(0);
01771 if (1==_protGenLevel) {
01772 protData=(RooDataSet*)_protDataset->reduce(fullProtVars);
01773 if (toyVerbose) {
01774 cout<<"Toy The prototype dataset for toy study:"<<endl;
01775 protData->Print("v");
01776
01777 cout<<endl;
01778 }
01779 }
01780
01781
01782 string fParamSStr;
01783 writeToStr(fullParams, fParamSStr);
01784
01785 Int_t nLoops=1;
01786 Int_t nExpPerLoop=toyNexp;
01787 Bool_t tooManyEvts=kFALSE;
01788 if (nExpPerLoop*toyNevt>1000000) tooManyEvts=kTRUE;
01789 if (tooManyEvts||(_protGenLevel>1)||(preToyRandParser.nArgs()>0)) {
01790 nLoops=toyNexp;
01791 nExpPerLoop=1;
01792 }
01793
01794 Bool_t firstToy(kTRUE);
01795
01796 for (Int_t expIdx=0; expIdx<nLoops; expIdx++) {
01797
01798
01799 TString toyFileName=Form(toyFilePrefix.Data(), expIdx);
01800 if (nExpPerLoop>1) toyFileName+=".%03d";
01801 toyFileName+=".text";
01802
01803 readFromStr(fullParams, fParamSStr);
01804 readFromStr(*fCloneSet, fParamSStr);
01805
01806 if ((preToyRandParser.nArgs()>0)&&_theToyParamGen) {
01807
01808 RooDataSet *theRandDataSet=_theToyParamGen->generate(paramRandSet, 1);
01809 RooArgSet *theRandSet=(RooArgSet*)theRandDataSet->get(0);
01810 stringstream theRandSetStr;
01811 theRandSet->writeToStream((ostream&)theRandSetStr, kFALSE);
01812 delete theRandDataSet;
01813 paramRandSet.readFromStream((istream&)theRandSetStr, kFALSE);
01814 cout<<"Toy Randomizer ArgSet:"<<endl;
01815 paramRandSet.Print("v");
01816 Int_t nParams=preToyRandParser.nArgs()/2;
01817 RooArgSet rParams, cParams;
01818 for (Int_t i=0; i<nParams; i++) {
01819 RooRealVar *theParam=dynamic_cast<RooRealVar*>
01820 (getAbsVar(preToyRandParser[2*i]));
01821 if (!theParam) {
01822 cout<<"Toy Can not find param "<<preToyRandParser[2*i]<<endl;
01823 continue;
01824 }
01825 rParams.add(*theParam);
01826
01827 theParam=(RooRealVar*)fCloneSet->find(theParam->GetName());
01828 if (!theParam) {
01829 cout<<"Toy Can not find cloned param "<<theParam->GetName()<<endl;
01830 exit(-1);
01831 }
01832 cParams.add(*theParam);
01833 cout<<"Toy "<<theParam->GetName()<<" = "<<theParam->getVal();
01834 theParam->setVal(getFormulaVal(preToyRandParser[2*i+1]));
01835 cout<<" --> "<<theParam->getVal()<<endl;
01836 }
01837
01838 string cTorStr;
01839 writeToStr(cParams, cTorStr);
01840 readFromStr(rParams, cTorStr);
01841 }
01842
01843 string randParamSStr;
01844 writeToStr(fullParams, randParamSStr);
01845
01846 TString genStr="";
01847 Double_t toyNevtGen=toyNevt;
01848 for (Int_t i=0; i<nEvtType; i++) {
01849 RooStringVar *evtTypeStr=(RooStringVar *)evtTypeStrList.at(i);
01850 rarStrParser toySrcParser=evtTypeStr->GetName();
01851 rarStrParser toySrcStrParser=evtTypeStr->getVal();
01852 if ("pdf"==toySrcParser[1]) continue;
01853
01854 Double_t evtTypeStrVal(0);
01855 while (toySrcStrParser.nArgs()>0) {
01856 evtTypeStrVal+=getFormulaVal(toySrcStrParser[0]);
01857 toySrcStrParser.Remove();
01858 }
01859
01860 Double_t nEvtWSrc(0);
01861 for (Int_t j=0; j<nCoeff; j++) {
01862 TString coefName=compCoeffList[j].GetName();
01863 nEvtWSrc+=((RooAbsReal*)fCloneSet->find(coefName))->getVal();
01864 }
01865 RooRealVar *paramVar=(RooRealVar*)fCloneSet->find(toySrcParser[0]);
01866 paramVar->setVal(paramVar->getVal()-evtTypeStrVal);
01867 Double_t nEvtWoSrc(0);
01868 for (Int_t j=0; j<nCoeff; j++) {
01869 TString coefName=compCoeffList[j].GetName();
01870 nEvtWoSrc+=((RooAbsReal*)fCloneSet->find(coefName))->getVal();
01871 }
01872 Double_t nEvtGen=nEvtWSrc-nEvtWoSrc;
01873 if (nEvtGen<0) {
01874 cout<<"Toy W A R N I N G !"<<endl
01875 <<" Can not generate negative event from src "<<toySrcParser[1]
01876 <<" for parameter "<<toySrcParser[0]<<endl;
01877 continue;
01878 }
01879 genStr+=Form(" \"%s\" %f ", toySrcParser[1].Data(), nEvtGen);
01880 toyNevtGen-=nEvtGen;
01881 }
01882 if (""!=genStr) {
01883 cout<<"Toy Using generation string: \""<<genStr<<"\" for embedding"<<endl;
01884 }
01885
01886 if (_protGenLevel>=2) {
01887 RooArgSet protDeps(fullProtVars);
01888
01889 protDeps.add(_compCat);
01890
01891 Int_t nEvt=_protDataset->numEntries();
01892 if (protData) delete protData;
01893 if (_protGenLevel>=3) protData=theProtGen->generate(protDeps, nEvt);
01894 else {
01895
01896 RooArgSet fCatSet(*((RooSimultaneous*)theProtGen)
01897 ->indexCat().getDependents(_protDataset));
01898 fCatSet.add(_compCat);
01899 RooDataSet *fCatData=theProtGen->generate(fCatSet, nEvt);
01900
01901 RooArgSet catSet(fCatSet);
01902
01903 catSet.remove(_protDataEVars, kFALSE, kTRUE);
01904 RooSuperCategory sCat("sCat", "sCat", catSet);
01905 sCat.attachDataSet(*fCatData);
01906 protData=new RooDataSet("protData", "protData", protDeps);
01907 for (Int_t i=0; i<nEvt; i++) {
01908 const RooArgSet *theCat=fCatData->get(i);
01909 TString dsName=sCat.getLabel();
01910 dsName.ReplaceAll("{", "");
01911 dsName.ReplaceAll("}", "");
01912
01913 RooDataSet *theData=(RooDataSet*)_protDatasets.FindObject(dsName);
01914 if (!theData) {
01915 if (dsName.Last(';')<0) theData=(RooDataSet*)_protDatasetsM.At(0);
01916 else {
01917 dsName.Remove(dsName.Last(';'), dsName.Length());
01918 theData=(RooDataSet*)_protDatasetsM.FindObject(dsName);
01919 }
01920 }
01921 if (!theData) theData=(RooDataSet*)_protDatasetsM.At(0);
01922 Int_t I=RooRandom::randomGenerator()->Integer(theData->numEntries());
01923 RooArgSet protSet(*theCat);
01924 protSet.add(*theData->get(I), kTRUE);
01925 protData->add(protSet);
01926 }
01927 delete fCatData;
01928 }
01929 cout<<"Toy protData from protGen:"<<endl;
01930 protData->Print("v");
01931
01932
01933
01934
01935 RooArgList protDepList(protDeps);
01936 for (Int_t i=0; i<protDepList.getSize(); i++) {
01937 RooCategory *theCat=(RooCategory *)&(protDepList[i]);
01938 if (theCat->ClassName()!=TString("RooCategory")) continue;
01939 Roo1DTable *theTable(0);
01940 if (theCat&&(theTable=protData->table(*theCat))) {
01941 theTable->Print();
01942
01943 TIterator* catTypeIter = theCat->typeIterator();
01944 RooCatType *theType(0);
01945 while(theType=(RooCatType*)catTypeIter->Next()) {
01946 TString typeName=theType->GetName();
01947 cout<<" "<<typeName<<"\t"<<theTable->getFrac(typeName)<<endl;
01948 }
01949 cout<<endl;
01950 delete catTypeIter;
01951 }
01952 }
01953 }
01954
01955 RooArgSet toyDeps(*theGen->getDependents(_protDataset));
01956 if (protData) toyDeps.remove(*theGen->getDependents(protData));
01957
01958 toyDeps.Print();
01959 RooMCStudy *theToy(0);
01960 if (toyFitMinosAS.getSize()<1) {
01961 theToy=new RooMCStudy
01962 (*theGen, toyDeps, FitModel(*_thePdf),
01963 Extended(extendedGen), ProtoData(*protData, kTRUE),
01964 ProjectedObservables(_conditionalObs),
01965 FitOptions(Save(kTRUE), Extended(kTRUE),
01966 Verbose(toyFitVerbose), Hesse(toyFitHesse),
01967 Minos(toyFitMinos), NumCPU(toyFitNumCPU)));
01968 } else {
01969 theToy=new RooMCStudy
01970 (*theGen, toyDeps, FitModel(*_thePdf),
01971 Extended(extendedGen), ProtoData(*protData, kTRUE),
01972 ProjectedObservables(_conditionalObs),
01973 FitOptions(Save(kTRUE), Extended(kTRUE),
01974 Verbose(toyFitVerbose), Hesse(toyFitHesse),
01975 Minos(toyFitMinos), Minos(toyFitMinosAS)), FitOptions(NumCPU(toyFitNumCPU)));
01976 }
01977 if (firstToy) {
01978 cout<<"Toy The toyDeps"<<endl;
01979 toyDeps.Print("v");
01980 cout<<"Toy The protVars"<<endl;
01981 _protDataVars.Print("v");
01982 }
01983
01984 Bool_t toyChkNegativePdf(kFALSE);
01985 if ("yes"==readConfStr("toyChkNegativePdf", "no", _runSec)) {
01986 toyChkNegativePdf=kTRUE;
01987 }
01988 RooDataSet *chkNPdfDS(0);
01989
01990 if ("no"!=readConfStr("toyGenerate", "yes", _runSec)) {
01991 if (firstToy) {
01992 cout<<endl<<"Toy RooMCStudy params for generating:"<<endl;
01993 theGen->getParameters(_protDataset)->Print("v");
01994 }
01995
01996 if (toyChkNegativePdf) {
01997 chkNPdfDS=_dummyPdf.generate(toyDeps, *protData, (Int_t)toyNevt);
01998 }
01999
02000 theToy->generate(nExpPerLoop, randInt(toyNevtGen), kTRUE);
02001 RooArgSet toyObs(toyDeps);
02002 if (protData) toyObs.add(*protData->get(),kTRUE);
02003 if (firstToy) {
02004 cout<<"Toy Full obs generated:"<<endl;
02005 toyObs.Print();
02006 }
02007 if (kTRUE) {
02008
02009 }
02010 if (""!=genStr) {
02011
02012 RooArgList etoyDeps(toyObs);
02013 etoyDeps.remove(_protDataEVars, kFALSE, kTRUE);
02014 if (firstToy) {
02015 cout<<"Toy etoyDeps"<<endl;
02016 etoyDeps.Print("v");
02017 cout<<"Toy protDataEVars"<<endl;
02018 _protDataEVars.Print("v");
02019 cout<<"Toy genStr "<<genStr<<endl;
02020 }
02021 generate(theToy, etoyDeps, genStr, genOpt, nExpPerLoop);
02022 }
02023 if (!toyFileName.BeginsWith("no")) {
02024 Int_t nSamples=nExpPerLoop;
02025 while (nSamples--) {
02026 TString thisToyFileName=Form(toyFileName.Data(),nSamples);
02027 ((RooDataSet*)theToy->genData(nSamples))->write(thisToyFileName);
02028
02029
02030 thisToyFileName.ReplaceAll(".text", ".root");
02031 #ifndef USENEWROOT
02032 TFile f(thisToyFileName, "recreate");
02033 ((RooDataSet*)theToy->genData(nSamples))->tree().Write();
02034 f.Close();
02035 #else
02036 RooDataSet *theSet = (RooDataSet*) theToy->genData(nSamples);
02037 saveAsRootFile(theSet, thisToyFileName, kTRUE);
02038 #endif
02039 }
02040 }
02041 if (firstToy) {
02042 RooDataSet *toySample=(RooDataSet*)theToy->genData(0);
02043 toySample=(RooDataSet*)toySample->Clone();
02044 toySample->SetName(_datasets->getDSName("toySample"));
02045 _datasets->getDatasetList()->Add(toySample);
02046 _datasets->ubStr(_datasets->getDSName("toySample"), "Unblinded");
02047 }
02048 }
02049
02050 if ("no"!=readConfStr("toyFit", "yes", _runSec)) {
02051 if (firstToy) {
02052 cout<<endl<<"Toy fitting coeffParamSet:"<<endl;
02053 coeffParamSet.Print("v");
02054 }
02055
02056 if (toyChkNegativePdf&&chkNPdfDS) {
02057
02058 attachDataSet(*chkNPdfDS);
02059 Int_t nChkEvt=chkNPdfDS->numEntries();
02060 for (Int_t iChkEvt=0; iChkEvt<nChkEvt; iChkEvt++) {
02061 RooArgSet *theEvt=(RooArgSet*)chkNPdfDS->get(iChkEvt);
02062 if (isNegativeValue()) {
02063 cout<<"Toy The PDF have negative value for current event"<<endl;
02064 theEvt->Print();
02065 theEvt->Print("v");
02066 exit(-1);
02067 }
02068 }
02069 }
02070 if ("no"!=readConfStr("toyGenerate", "yes", _runSec)) {
02071
02072 TList genSamples;
02073 Int_t nSamples=nExpPerLoop;
02074 while (nSamples--) genSamples.Add(theToy->genData(nSamples)->Clone());
02075 theToy->fit(nExpPerLoop, genSamples);
02076 } else {
02077 theToy->fit(nExpPerLoop, toyFileName);
02078 }
02079
02080 if (""!=genStr) {
02081 cout<<"Toy Re-calculate pulls for embedded toys"<<endl;
02082
02083 readFromStr(fullParams, randParamSStr);
02084 RooDataSet &fitParData=(RooDataSet &)theToy->fitParDataSet();
02085 const RooArgSet *fitParDataSet=fitParData.get();
02086 TIterator* iter = coeffParamSet.createIterator();
02087 RooRealVar *par(0);
02088 while(par=(RooRealVar*)iter->Next()) {
02089 TString pullName=Form("%spull", par->GetName());
02090 RooPullVar *origPull=(RooPullVar*)fitParDataSet->find(pullName);
02091 if (!origPull) continue;
02092 RooAbsReal* tPar=(RooAbsReal*)par->Clone("truth");
02093 RooPullVar pull(pullName+"_embd", origPull->GetTitle(),*par,*tPar);
02094 fitParData.addColumn(pull);
02095 delete tPar;
02096 }
02097 delete iter;
02098
02099 rarStrParser genStrParser=genStr;
02100 Int_t nSubset=genStrParser.nArgs()/2;
02101 for(Int_t i=0; i<nSubset; i++) {
02102 TString genSrcName=genStrParser[0];
02103 genStrParser.Remove();
02104 Double_t nEvtGen=atof(genStrParser[0]);
02105 genStrParser.Remove();
02106 RooRealVar embdEvt("embdEvt_"+genSrcName, "embd Evt", nEvtGen);
02107 fitParData.addColumn(embdEvt);
02108 }
02109 }
02110 {
02111 RooDataSet &fitParData=(RooDataSet &)theToy->fitParDataSet();
02112 RooRealVar ID_COL0("ID_COL0", "Toy ID", _toyID);
02113 fitParData.addColumn(ID_COL0);
02114 RooRealVar ID_COL1("ID_COL1", "Toy Loop ID", expIdx);
02115 fitParData.addColumn(ID_COL1);
02116 RooRealVar ID_COL2("ID_COL2", "Toy Loop #", 0);
02117 fitParData.addColumn(ID_COL2);
02118
02119 const RooFitResult *fr=theToy->fitResult(0);
02120 Int_t nFloat=fr->floatParsFinal().getSize();
02121 for (Int_t iF=1; iF<=nFloat; iF++) {
02122 RooRealVar CorrG(Form("CorrG_%i",iF), "Global Corr",0);
02123 fitParData.addColumn(CorrG);
02124 for(Int_t jF=1; jF<iF; jF++) {
02125 RooRealVar Corr(Form("Corr_%i_%i", iF, jF), "Corr",0);
02126 fitParData.addColumn(Corr);
02127 }
02128 }
02129
02130 RooRealVar covQual("covQual", "covQual", 0);
02131 fitParData.addColumn(covQual);
02132 RooRealVar numInvalidNLL("numInvalidNLL", "numInvalidNLL", 0);
02133 fitParData.addColumn(numInvalidNLL);
02134 RooRealVar edm("edm", "edm", 0);
02135 fitParData.addColumn(edm);
02136
02137
02138 RooRealVar GOFChisq("GOFChisq", "GOFChisq", 0);
02139 fitParData.addColumn(GOFChisq);
02140 }
02141 if (!toyResults) {
02142 toyResults=new
02143 RooDataSet("toyResults","toyResults",*theToy->fitParDataSet().get());
02144 }
02145
02146 Int_t ii=0;
02147 for (Int_t i=0; i<nExpPerLoop; i++) {
02148 RooFitResult *fr=(RooFitResult*)theToy->fitResult(i);
02149 if (fr->status()) {
02150 cout<<"Toy Fit status for experiment #"<<expIdx<<"-"<<nExpPerLoop-i
02151 <<": "<<fr->status()<<endl;
02152 continue;
02153 }
02154 RooArgSet *fitResultSet=(RooArgSet *)theToy->fitParams(ii);
02155 ((RooRealVar*)fitResultSet->find("ID_COL2"))->setVal(nExpPerLoop-i);
02156
02157 Int_t nFloat=fr->floatParsFinal().getSize();
02158 for (Int_t iF=0; iF<nFloat; iF++) {
02159 RooAbsArg &iArg=fr->floatParsFinal().operator[](iF);
02160 ((RooRealVar*)fitResultSet->find(Form("CorrG_%i",iF+1)))->
02161 setVal(fr->globalCorr(iArg));
02162 for(Int_t jF=0; jF<iF; jF++) {
02163 RooAbsArg &jArg=fr->floatParsFinal().operator[](jF);
02164 ((RooRealVar*)fitResultSet->find(Form("Corr_%i_%i", iF+1,jF+1)))->
02165 setVal(fr->correlation(iArg, jArg));
02166 }
02167 }
02168
02169 ((RooRealVar*)fitResultSet->find("covQual"))->setVal(fr->covQual());
02170 ((RooRealVar*)fitResultSet->find("numInvalidNLL"))
02171 ->setVal(fr->numInvalidNLL());
02172 ((RooRealVar*)fitResultSet->find("edm"))->setVal(fr->edm());
02173
02174
02175 Double_t gofChisq(0);
02176 TString postMLGOFChisq=readConfStrCnA("postMLGOFChisq", "no");
02177 if (!postMLGOFChisq.BeginsWith("no"))
02178 gofChisq=doGOFChisq((RooDataSet*)theToy->genData(i), cout);
02179 ((RooRealVar*)fitResultSet->find("GOFChisq"))->setVal(gofChisq);
02180
02181 toyResults->add(*fitResultSet);
02182 ii++;
02183 }
02184 }
02185 delete theToy;
02186 firstToy=kFALSE;
02187 }
02188
02189
02190 return toyResults;
02191 }
02192
02199 void rarMLFitter::getCompCatDS(TList*ds,RooDataSet *iData,RooCategory *compCat)
02200 {
02201 RooArgSet catSet(*((RooSimultaneous*)_theProtGen)
02202 ->indexCat().getDependents(iData));
02203 catSet.remove(_compCat, kFALSE, kTRUE);
02204 catSet.remove(_protDataEVars, kFALSE, kTRUE);
02205 if (catSet.getSize()<=0) {
02206 RooDataSet *theData=(RooDataSet*)iData->reduce("1");
02207 theData->SetName("noCompCatDS");
02208 if (compCat) {
02209 theData->addColumn(*compCat);
02210 theData->SetName(compCat->getLabel());
02211 }
02212 ds->Add(theData);
02213 return;
02214 }
02215
02216 RooSuperCategory sCat("sCat", "sCat", catSet);
02217 sCat.attachDataSet(*iData);
02218 Int_t nEvt=iData->numEntries();
02219 for (Int_t i=0; i<nEvt; i++) {
02220 RooArgSet *theEvt=(RooArgSet*)iData->get(i);
02221 TString dsName=sCat.getLabel();
02222 if (compCat) dsName=Form("%s;%s", dsName.Data(), compCat->getLabel());
02223 dsName.ReplaceAll("{", "");
02224 dsName.ReplaceAll("}", "");
02225 RooDataSet *theData=(RooDataSet*)ds->FindObject(dsName);
02226 if (!theData) {
02227 theData=new RooDataSet(dsName, dsName, *theEvt);
02228 if (compCat) theData->addColumn(*compCat);
02229 ds->Add(theData);
02230 }
02231 theData->add(*theEvt);
02232 }
02233 }
02234
02241 Int_t rarMLFitter::randInt(Double_t iNumber)
02242 {
02243 Int_t retVal=(Int_t) iNumber;
02244 if (RooRandom::randomGenerator()->Uniform()+retVal<iNumber) retVal++;
02245 return retVal;
02246 }
02247
02249 void rarMLFitter::generate(RooMCStudy *theToy, RooAbsPdf *theGen,
02250 const TString genOpt, const Int_t toyNexp)
02251 {
02252 cout<<endl<<" In rarMLFitter generate (two-step) for "<<GetName()<<endl;
02253
02254
02255 Int_t nSamples=toyNexp;
02256 while(nSamples--) {
02257 cout<<endl;
02258 RooDataSet *genSample=(RooDataSet *)theToy->genData(nSamples);
02259
02260
02261
02262 RooArgSet subSet(*genSample->get());
02263 subSet.remove(*subSet.selectByName("hzCat,hcCat"),
02264 kTRUE, kTRUE);
02265
02266 RooDataSet *subData=(RooDataSet*)genSample->reduce(subSet);
02267
02268 RooArgSet *addOns=getAddOnCols();
02269 addOns=(RooArgSet*)(addOns->selectByName("hzCat,hcCat"));
02270 subData->addColumns(*addOns);
02271
02272
02273 genSample->reset();
02274 genSample->append(*subData);
02275
02276 delete subData;
02277 }
02278 if (genOpt.Contains("safasse")) {theGen->Print();}
02279 }
02280
02294 void rarMLFitter::generate(RooMCStudy *theToy,
02295 const RooArgList& etoyDeps,
02296 const TString genStr, const TString genOpt,
02297 const Int_t toyNexp)
02298 {
02299 cout<<endl<<" In rarMLFitter generate (embed) for "<<GetName()<<endl;
02300
02301
02302 Int_t nSamples=toyNexp;
02303 while(nSamples--) {
02304 cout<<endl;
02305 RooDataSet *genSample=(RooDataSet *)theToy->genData(nSamples);
02306
02307 rarStrParser genStrParser=genStr;
02308 Int_t nSubset=genStrParser.nArgs()/2;
02309 for(Int_t i=0; i<nSubset; i++) {
02310 TString genSrcName=genStrParser[0];
02311 genStrParser.Remove();
02312 Double_t nEvtGen=atof(genStrParser[0]);
02313 genStrParser.Remove();
02314 RooDataSet *embedSample=generate(etoyDeps, genSrcName, nEvtGen, genOpt);
02315 if (embedSample) {
02316 genSample->append(*embedSample);
02317 delete embedSample;
02318 }
02319 }
02320 }
02321 }
02322
02338 RooDataSet *rarMLFitter::generate(const RooArgList& dependents,
02339 const TString genSrcName, Double_t nEvtGen,
02340 const TString genOpt)
02341 {
02342
02343 RooDataSet *theSample(0);
02344 RooArgSet fullDeps(dependents);
02345 RooArgSet randObs;
02346 RooAbsPdf *theObsGen(0);
02347
02348 RooStringVar *srcRandStr=(RooStringVar *)_embdObsRandSet.find(genSrcName);
02349 if (srcRandStr) {
02350 rarStrParser srcRandStrParser=srcRandStr->getVal();
02351
02352 while (srcRandStrParser.nArgs()>0) {
02353 TString obsName=srcRandStrParser[0];
02354 srcRandStrParser.Remove();
02355 if (_fObsSet.find(obsName)) randObs.add(*_fObsSet.find(obsName));
02356 else theObsGen=(RooAbsPdf*)_embdObsGens.find("the_"+obsName);
02357 }
02358 if (theObsGen) {
02359 cout<<" Using "<<theObsGen->GetName()<<" as embd randomizer"<<endl;
02360 theObsGen->Print();
02361
02362 } else {
02363 cout<<" Can not find randomizer defined in \""<<srcRandStr->getVal()
02364 <<"\""<<endl;
02365 }
02366 if (randObs.getSize()>0) {
02367 cout<<" The following params will be randomized:"<<endl;
02368 randObs.Print("v");
02369 } else {
02370 cout<<" Can not find params to randomize in \""<<srcRandStr->getVal()
02371 <<"\""<<endl;
02372 }
02373 }
02374 if (theObsGen) fullDeps.add(*theObsGen->getDependents(_protDataset), kTRUE);
02375
02376 fullDeps.remove(_compCat, kFALSE, kTRUE);
02377
02378 RooDataSet *genSrcData=_datasets->getData(genSrcName);
02379 if (!genSrcData) {
02380 cout<<"toy src data \""<<genSrcName<<"\" does not exist"<<endl;
02381 exit(-1);
02382 }
02383
02384 if (genSrcData->numEntries()<=0) {
02385 cout<<"toy src data \""<<genSrcName<<"\" does not contain any entries"
02386 <<endl;
02387 exit(-1);
02388 }
02389
02390 if (genOpt.Contains("e")) {
02391 cout<<"Extended generating nEvt "<<nEvtGen;
02392 nEvtGen = RooRandom::randomGenerator()->Poisson(nEvtGen);
02393 cout<<" -> "<<nEvtGen<<endl;
02394 }
02395 nEvtGen=randInt(nEvtGen);
02396 cout<<"Generate "<<nEvtGen<<" events from dataset "
02397 <<genSrcData->GetName()<<endl;
02398 if (nEvtGen<.1) {
02399 return theSample;
02400 }
02401
02402 rarStrParser embdUnCorrParser=readConfStr("toyEmbdUnCorrelate","",_runSec);
02403 for (Int_t i=0; i<=embdUnCorrParser.nArgs(); i++) {
02404 RooDataSet *subSample(0);
02405 if ((i==embdUnCorrParser.nArgs())&&(fullDeps.getSize()>0)) {
02406
02407
02408 subSample=new RooDataSet("subSample", "subSample", fullDeps);
02409 } else {
02410 RooArgSet subSampleSet;
02411 rarStrParser groupObsParser=embdUnCorrParser[i];
02412 for (Int_t j=0; j<groupObsParser.nArgs(); j++) {
02413 RooAbsArg *theObs=fullDeps.find(groupObsParser[j]);
02414 if (theObs) {
02415 subSampleSet.add(*theObs);
02416 fullDeps.remove(*theObs);
02417 }
02418 }
02419 if (subSampleSet.getSize()<1) continue;
02420 cout<<" Sub obs:"<<endl;
02421 subSampleSet.Print();
02422 subSample=new RooDataSet("subSample", "subSample", subSampleSet);
02423 }
02424 if(!subSample) continue;
02425 for (Int_t j=0; j<nEvtGen; j++) {
02426 Int_t I=RooRandom::randomGenerator()->Integer(genSrcData->numEntries());
02427 const RooArgSet* row=genSrcData->get(I);
02428 if(!row) {
02429 cout<<"skipping null row entry I="<<I<<endl;
02430 continue;
02431 }
02432 subSample->add(*row);
02433 }
02434 if (!theSample)
02435 theSample=subSample;
02436 else {
02437 theSample->merge(subSample);
02438 delete subSample;
02439 }
02440 }
02441
02442
02443 if (_protDataEVars.getSize()>0) {
02444 RooDataSet *protSample=
02445 new RooDataSet("protEDataSet", "prot embedded Dataset", _protDataEVars);
02446 Int_t nEvt=theSample->numEntries();
02447 for (Int_t i=0; i<nEvt; i++) {
02448 Int_t Idx=RooRandom::randomGenerator()->
02449 Integer(_protDataset->numEntries());
02450 const RooArgSet* row=_protDataset->get(Idx);
02451 protSample->add(*row);
02452 }
02453 cout<<"Merge "<<nEvt<<" events from prototype dataset "
02454 <<_protDataset->GetName()<<endl;
02455 theSample->merge(protSample);
02456 delete protSample;
02457 }
02458 {
02459 _compCat.setIndex(_compCat.numTypes()-1);
02460 theSample->addColumn(_compCat);
02461 }
02462
02463 if (theObsGen&&randObs.getSize()>0) {
02464 RooArgSet embProtSet;
02465 embProtSet.add(*theSample->get(0));
02466 embProtSet.remove(randObs, kFALSE, kTRUE);
02467 RooDataSet *protDS=(RooDataSet*)theSample->reduce(embProtSet);
02468 RooDataSet *oldembedSample=theSample;
02469 theSample=theObsGen->generate(randObs,*protDS, (Int_t)nEvtGen);
02470
02471
02472
02473
02474
02475
02476 delete oldembedSample;
02477 delete protDS;
02478 }
02479
02480 return theSample;
02481 }
02482
02487 RooAbsPdf *rarMLFitter::getProtGen()
02488 {
02489 if (_theProtGen) return _theProtGen;
02490
02491 rarCompBase::getProtGen();
02492
02493 if (getControlBit("noSimFit")) {
02494 _theProtGen=(RooAbsPdf*)_subProtGenPdfs.at(0);
02495 } else {
02496 RooSimPdfBuilder *simBuilder=new RooSimPdfBuilder(_subProtGenPdfs);
02497 RooArgSet *simConfig=simBuilder->createProtoBuildConfig();
02498
02499 ((RooStringVar*)simConfig->find("splitCats"))->
02500 setVal(((RooStringVar*)_simConfig->find("splitCats"))->getVal());
02501 TString modelStr=((RooStringVar*)_simConfig->find("physModels"))->getVal();
02502 for (Int_t i=0; i<_subProtGenPdfs.getSize(); i++) {
02503 TString genName = _subProtGenPdfs[i].GetName();
02504 TString pdfName =_subPdfs[i].GetName();
02505
02506 TString configStr=((RooStringVar*)_simConfig->find(pdfName))->getVal();
02507 ((RooStringVar*)simConfig->find(genName))->setVal(configStr);
02508 modelStr.ReplaceAll(pdfName, genName);
02509 }
02510 ((RooStringVar*)simConfig->find("physModels"))->setVal(modelStr);
02511 cout<<"simPdfBuilder configs for protGen:"<<endl;
02512 simConfig->Print("v");
02513 simBuilder->addSpecializations(_simBuilder->splitLeafList());
02514 _theProtGen=(RooAbsPdf*)
02515 simBuilder->buildPdf(*simConfig, _protDataset, 0, kTRUE);
02516
02517 RooAbsCategoryLValue *simCats=(RooAbsCategoryLValue *)
02518 (&((RooSimultaneous*)_theProtGen)->indexCat());
02519 Int_t nCats=simCats->numTypes();
02520 for(Int_t i=0; i<nCats; i++) {
02521 TString catName=((RooCatType*)simCats->lookupType(i))->GetName();
02522 if (!((RooSimultaneous*)_theProtGen)->getPdf(catName.Data())) {
02523 ((RooSimultaneous*)_theProtGen)->addPdf(_dummyExtPdf, catName);
02524 }
02525 }
02526
02527 _theProtGen->SetName(Form("simProtGen_%s", GetName()));
02528 _theProtGen->Print();
02529
02530 }
02531
02532 _theProtGen=compGen(_theProtGen, _subProtGenPdfs, _compCat);
02533
02534 return _theProtGen;
02535 }
02536
02541 RooAbsPdf *rarMLFitter::getGenerator()
02542 {
02543 if (_theGen) return _theGen;
02544 _theGen=_thePdf;
02545 if (getControlBit("SimFit")) {
02546 RooAbsCategoryLValue *simCats=(RooAbsCategoryLValue *)
02547 (&((RooSimultaneous*)_thePdf)->indexCat());
02548 Int_t nCats=simCats->numTypes();
02549
02550 _theGen=new RooSimultaneous
02551 (Form("simGen_%s",GetName()),"simGen",*simCats);
02552 for(Int_t i=0; i<nCats; i++) {
02553 TString catName=((RooCatType*)simCats->lookupType(i))->GetName();
02554 RooAbsPdf *catPdf=((RooSimultaneous*)_thePdf)->getPdf(catName.Data());
02555 if (!catPdf) catPdf=&_dummyExtPdf;
02556 ((RooSimultaneous*)_theGen)->addPdf(*catPdf, catName);
02557 }
02558 }
02559
02560 if (_protGenLevel>1) _theGen=compGen(_theGen, _subPdfs, _compCat);
02561
02562 cout<<"Generator created from fitter:"<<endl;
02563 _theGen->Print();
02564 return _theGen;
02565 }
02566
02575 RooAbsArg *rarMLFitter::findSimed(RooArgSet &simSet,
02576 TString argName, TString catName,
02577 RooAbsPdf *srcPdf)
02578 {
02579 RooAbsArg *theComp(0);
02580
02581 theComp=simSet.find(argName+"_"+catName);
02582 if (theComp) return theComp;
02583
02584 if (catName.Contains(";")) {
02585
02586 TString phyCat=catName(1, catName.First(";")-1);
02587 catName.Replace(1, phyCat.Length()+1,"");
02588 catName.Replace(catName.Length()-1, 0, ";"+phyCat);
02589 theComp=simSet.find(argName+"_"+catName);
02590 if (theComp) return theComp;
02591 }
02592
02593 theComp=simSet.find(argName);
02594 if (""==getPhysCat()||" "==getPhysCat()) return theComp;
02595 if (!srcPdf) return theComp;
02596 RooSimultaneous *sPdf=dynamic_cast<RooSimultaneous*>(srcPdf);
02597 if (!sPdf) return theComp;
02598
02599 theComp=0;
02600 RooAbsPdf *mlPdf=sPdf->getPdf(catName);
02601 if (!mlPdf) return theComp;
02602 RooArgSet *mlSet=mlPdf->getComponents();
02603 theComp=mlSet->find(argName);
02604 return theComp;
02605 }
02606
02615 RooSimultaneous *rarMLFitter::compGen(RooAbsPdf *gen, RooArgList subGens,
02616 RooCategory &compCat)
02617 {
02618 RooSimultaneous *theCompGen(0);
02619 if (getControlBit("noSimFit")) {
02620 Int_t nCompTypes=compCat.numTypes();
02621 theCompGen=new RooSimultaneous
02622 (Form("simComp_%s", gen->GetName()), "comp-cated gen", compCat);
02623 RooAddPdf &theMode=(RooAddPdf&)subGens[0];
02624 RooArgList pdfList(theMode.pdfList());
02625 RooArgList coefList(theMode.coefList());
02626 Int_t nComps=pdfList.getSize();
02627 for (Int_t i=0; i<nComps; i++) {
02628 RooAddPdf *thePdf=new RooAddPdf
02629 (Form("comp_%s_%s", theMode.GetName(), pdfList[i].GetName()),
02630 Form("comp %s %s", theMode.GetName(), pdfList[i].GetName()),
02631 pdfList[i], coefList[i]);
02632 theCompGen->addPdf(*thePdf, Form("%02d", i));
02633 }
02634
02635 for (Int_t i=nComps; i<nCompTypes; i++)
02636 theCompGen->addPdf(_dummyExtPdf, Form("%02d", i));
02637
02638 return theCompGen;
02639 }
02640
02641 RooAbsCategoryLValue *simCats=(RooAbsCategoryLValue *)
02642 (&((RooSimultaneous*)gen)->indexCat());
02643 RooArgSet *genSimSet=gen->getComponents();
02644
02645 RooArgSet myInputCats(compCat);
02646
02647 RooSuperCategory *simSuperCats=dynamic_cast<RooSuperCategory*>(simCats);
02648 if (simSuperCats) myInputCats.add(simSuperCats->inputCatList());
02649 else myInputCats.add(*simCats);
02650
02651 RooSuperCategory *compGenCat=new RooSuperCategory
02652 ("compGenCat", "compGenCat", myInputCats);
02653
02654 theCompGen=new RooSimultaneous
02655 (Form("simComp_%s", gen->GetName()), "comp-cated gen", *compGenCat);
02656 Int_t nCats=compGenCat->numTypes();
02657 for (Int_t i=0; i<nCats; i++) {
02658 TString catName=((RooCatType*)compGenCat->lookupType(i))->GetName();
02659 TString compLabel=catName(1,catName.First(";")-1);
02660 Int_t compIdx=atoi(compLabel);
02661
02662 TString gCatName=catName
02663 (catName.First(";")+1, catName.Last('}')-catName.First(";")-1);
02664 if (simSuperCats) gCatName="{"+gCatName+"}";
02665
02666 RooAbsPdf *thisCatGen(0);
02667 for (Int_t j=0; j<_nComp; j++) {
02668 RooAddPdf &theMode=(RooAddPdf&)subGens[j];
02669 if (compIdx>=theMode.pdfList().getSize()) continue;
02670
02671 TString thePdfName=theMode.pdfList()[compIdx].GetName();
02672 TString theCoefName=theMode.coefList()[compIdx].GetName();
02673 RooAbsPdf *thisPdf=(RooAbsPdf*)findSimed(*genSimSet,thePdfName,gCatName,
02674 gen);
02675 RooAbsReal *thisCoef=(RooAbsReal*)
02676 findSimed(*genSimSet, theCoefName, gCatName, gen);
02677 if (thisPdf&&thisCoef) {
02678 thisCatGen=new RooAddPdf
02679 ("comp_"+catName+"_"+thePdfName+"_"+theCoefName,
02680 "comp "+catName+" "+thePdfName+" "+theCoefName, *thisPdf,*thisCoef);
02681 break;
02682 }
02683 }
02684 if (thisCatGen) theCompGen->addPdf(*thisCatGen, catName);
02685 else theCompGen->addPdf(_dummyExtPdf, catName);
02686 }
02687
02688 theCompGen->Print();
02689
02690 RooArgSet nullDS;
02691 cout << "Expected number of events generated : "
02692 << theCompGen->expectedEvents(&nullDS)<<endl;
02693 for (Int_t i=0; i<nCats; i++) {
02694 TString catName=((RooCatType*)compGenCat->lookupType(i))->GetName();
02695 compGenCat->setLabel(catName);
02696 cout<<"\t"<<catName<<" "
02697 <<theCompGen->expectedEvents(&nullDS)<<endl;
02698 }
02699
02700 return theCompGen;
02701 }
02702
02707 RooAbsPdf *rarMLFitter::getExtCompPdf(RooAbsPdf *thePdf, RooAbsReal *theCoef)
02708 {
02709 RooAbsPdf *theComp(0);
02710 TString thePdfName=thePdf->GetName();
02711 TString theCoefName=theCoef->GetName();
02712 TString theCompName="ExtComp_"+thePdfName+"_"+theCoefName;
02713 if (theComp=(RooAbsPdf*)_extCompPdfs.find(theCompName)) return theComp;
02714
02715 if (!getControlBit("SimFit")) {
02716 theComp=new RooAddPdf
02717 (theCompName, theCompName, RooArgList(*thePdf), RooArgList(*theCoef));
02718 _extCompPdfs.add(*theComp);
02719 return theComp;
02720 }
02721
02722 RooAbsCategoryLValue *theCat=(RooAbsCategoryLValue *)
02723 (&((RooSimultaneous*)_thePdf)->indexCat());
02724 RooArgSet *theSimSet=_thePdf->getComponents();
02725 theComp=new RooSimultaneous(theCompName, theCompName, *theCat);
02726 _extCompPdfs.add(*theComp);
02727 TIterator* catTypeIter = theCat->typeIterator();
02728 RooCatType *theType(0);
02729 while(theType=(RooCatType*)catTypeIter->Next()) {
02730 TString theTypeName=theType->GetName();
02731 RooAbsPdf *thisPdf=(RooAbsPdf*)
02732 findSimed(*theSimSet, thePdfName, theTypeName);
02733 RooAbsReal *thisCoef=(RooAbsReal*)
02734 findSimed(*theSimSet, theCoefName, theTypeName);
02735 RooAbsPdf *thisCatPdf=&_dummyExtPdf;
02736 if (thisPdf&&thisCoef) {
02737 thisCatPdf=new RooAddPdf
02738 ("extComp_"+theTypeName+"_"+thePdfName+"_"+theCoefName,
02739 "extComp "+theTypeName+" "+thePdfName+" "+theCoefName,
02740 *thisPdf, *thisCoef);
02741 }
02742
02743
02744 ((RooSimultaneous *)theComp)->addPdf(*thisCatPdf, theTypeName);
02745 }
02746
02747 delete catTypeIter;
02748 return theComp;
02749 }
02750
02754 void rarMLFitter::getSnB()
02755 {
02756 if (_theSPdf) return;
02757
02758 rarStrParser compStrParser=readConfStr("projComps", "", _runSec);
02759 if (compStrParser.nArgs()<=0) {
02760 cout<<"no projection component defined with config \"projComps\""
02761 <<" in section \""<<_runSec<<"\""<<endl;
02762 exit(-1);
02763 }
02764
02765 RooArgList sList, bList;
02766 Int_t nComp=1;
02767 if (getControlBit("SimFit")) nComp=_nComp;
02768 for (Int_t i=0; i<nComp; i++) {
02769 rarMLPdf *mlPdf=(rarMLPdf*)_pdfList.At(i);
02770 TList *pdfList=mlPdf->getPdfList();
02771 Int_t pdfSize=pdfList->GetSize();
02772 RooArgList coefList(mlPdf->getCoeffList());
02773 for (Int_t j=0; j<pdfSize; j++) {
02774 rarBasePdf *thePdf=(rarBasePdf*)pdfList->At(j);
02775 RooAbsReal *theCoef=(RooAbsReal*)coefList.at(j);
02776 _fCompList.Add(thePdf);
02777 _fCoefList.add(*theCoef);
02778 RooAbsPdf *theExtCompPdf=getExtCompPdf(thePdf->getPdf(), theCoef);
02779
02780 if (compStrParser.Have(thePdf->GetName())) {
02781 sList.add(*theExtCompPdf);
02782 _sCompList.Add(thePdf);
02783 _sCoefList.add(*theCoef);
02784 } else bList.add(*theExtCompPdf);
02785 }
02786 }
02787 if ((sList.getSize()<1)||(bList.getSize()<1)) {
02788 cout<<"Can not find any comp for projPlot"<<endl;
02789 exit(-1);
02790 }
02791
02792 _theSPdf=new RooAddPdf("S_Pdf", "S Pdf", sList);
02793 _theBPdf=new RooAddPdf("B_Pdf", "B Pdf", bList);
02794
02795 return;
02796 }
02797
02804 RooDataSet *rarMLFitter::getLLRDataset(RooDataSet *theData, RooFormulaVar *LLR)
02805 {
02806 RooDataSet *retDS=new
02807 RooDataSet(*theData,Form("%s_wLLR", theData->GetName()));
02808
02809
02810
02811 setColLimits(retDS);
02812
02813
02814
02815
02816
02817 retDS->addColumn(*LLR);
02818
02819 setColLimits(retDS, kFALSE);
02820
02821
02822
02823
02824
02825
02826 return retDS;
02827 }
02828
02835 void rarMLFitter::doLLRPlot(RooDataSet *projData, TList &plotList)
02836 {
02837 cout<<endl<<" In rarMLFitter doLLRPlot for "<<GetName()<<endl;
02838 TString llrStr=readConfStr("projLLRPlots", "yes", _runSec);
02839 if ("no"==llrStr) return;
02840
02841 RooArgSet projVarSet;
02842 RooArgSet depVarSet(_fObsSet);
02843 depVarSet.remove(_conditionalObs,kFALSE,kTRUE);
02844
02845 const RooAbsReal *theSProj=_theSPdf->createPlotProjection(depVarSet, projVarSet);
02846 const RooAbsReal *theBProj=_theBPdf->createPlotProjection(depVarSet, projVarSet);
02847 RooArgList projPdfList(*theSProj, *theBProj, "Projected Pdf List");
02848 TString formulaStr="@0/(@0+@1)";
02849 TString funcName="LLR";
02850 RooFormulaVar *lRatioFunc=new
02851 RooFormulaVar(funcName, "LRatio Func", formulaStr, projPdfList);
02852
02853 Int_t nBins=atoi(readConfStr("plotBins_LLR", "100", _runSec));
02854
02855 TH1F *LLRP0=new TH1F(Form("LLR_%s", projData->GetName()),"LLR ds",nBins,0,1);
02856 projData=getLLRDataset(projData, lRatioFunc);
02857 for (Int_t idxEvt=0; idxEvt<projData->numEntries(); idxEvt++) {
02858 RooArgSet *dEntry=(RooArgSet*)projData->get(idxEvt);
02859 LLRP0->Fill(((RooRealVar*)dEntry->find(funcName))->getVal(),
02860 projData->weight());
02861 }
02862 LLRP0->Print("v");
02863 plotList.Add(LLRP0);
02864 RooDataSet *protData(0);
02866 if (_protDataVars.getSize()>0)
02867 protData=(RooDataSet*)projData->reduce(_protDataVars);
02868 RooArgSet theDeps(*_thePdf->getDependents(projData));
02869 if (protData) theDeps.remove(*_thePdf->getDependents(protData));
02870
02871
02872 TString projLLRScale=
02873 readConfStr(Form("projLLRScale_%s", GetName()), "-1", _runSec);
02874 if ("-1"==projLLRScale)
02875 projLLRScale=readConfStr("projLLRScale", "10", _runSec);
02876 Double_t scale=atof(projLLRScale);
02877 Int_t nEvt=(Int_t)(scale*projData->numEntries());
02878 if (nEvt>0) {
02879 getLLRPlot(plotList, GetName(), _thePdf, nEvt,
02880 theDeps, protData, nBins, lRatioFunc);
02881 } else {
02882 cout<<" Do not generate sample for LLR plot of "<<GetName()<<endl;
02883 }
02884
02885 TString llrCompStr="";
02886 if ("yes"!=llrStr) llrCompStr=llrStr;
02887 else for (Int_t i=0; i<_fCompList.GetSize(); i++)
02888 llrCompStr+=Form(" \"%s\"", _fCompList.At(i)->GetName());
02889 rarStrParser llrStrParser=llrCompStr;
02890 while (llrStrParser.nArgs()>0) {
02891 TString compName=llrStrParser[0];
02892 llrStrParser.Remove();
02893 rarBasePdf *compPdf=(rarBasePdf*)_fCompList.FindObject(compName);
02894 if (!compPdf) continue;
02895 Int_t compIdx=_fCompList.IndexOf(compPdf);
02896 RooAbsPdf *thePdf=compPdf->getSimPdf();
02897 if (!thePdf) thePdf=compPdf->getPdf();
02898
02899 TString projLLRScale=
02900 readConfStr(Form("projLLRScale_%s", compPdf->GetName()), "-1", _runSec);
02901 if ("-1"==projLLRScale)
02902 projLLRScale=readConfStr("projLLRScale", "10", _runSec);
02903 Double_t scale=atof(projLLRScale);
02904
02905 Int_t nEvt=(Int_t)(scale*((RooAbsReal&)_fCoefList[compIdx]).getVal());
02906 if (nEvt<=0) {
02907 cout<<" Negative yield for component "<<compPdf->GetName()
02908 <<", Ignored!!!"<<endl;
02909 continue;
02910 }
02911 getLLRPlot(plotList, compPdf->GetName(), thePdf, nEvt,
02912 theDeps, protData, nBins, lRatioFunc);
02913 }
02914
02915 delete projData;
02916
02917 return;
02918 }
02919
02931 void rarMLFitter::getLLRPlot(TList &plotList, TString plotName,
02932 RooAbsPdf *thePdf, Int_t nEvt,
02933 RooArgSet theDeps, RooDataSet *protData,
02934 Int_t nBins, RooAbsReal *LLRFunc)
02935 {
02936 TString histName="LLR_"+plotName;
02937 if (plotList.FindObject(histName)) return;
02938 TH1F *theHist(0);
02939
02940 RooDataSet *theData(0);
02941 if (protData)
02942 theData=thePdf->generate(theDeps,*protData,nEvt,kFALSE,kTRUE);
02943 else theData=thePdf->generate(theDeps, nEvt);
02944 theData->addColumn(*LLRFunc);
02945 theHist=new TH1F(histName, "LLR hist", nBins,0,1);
02946 for (Int_t idxEvt=0; idxEvt<theData->numEntries(); idxEvt++) {
02947 RooArgSet *dEntry=(RooArgSet*)theData->get(idxEvt);
02948 theHist->Fill(((RooRealVar*)dEntry->find(LLRFunc->GetName()))->getVal(),
02949 theData->weight());
02950 }
02951 theHist->Print("v");
02952 plotList.Add(theHist);
02953 delete theData;
02954 return;
02955 }
02956
02957
02969 RooPlot *rarMLFitter::doProjPlot(RooDataSet *projData, RooRealVar *theVar,
02970 TList &plotList)
02971 {
02972 cout<<endl<<" In rarMLFitter doProjPlot of "<<theVar->GetName()
02973 <<" for "<<GetName()<<endl;
02974 RooPlot *frame(0);
02975
02976 TString varName=theVar->GetName();
02977
02978 TString ppdStr=readConfStrCnA("projPlotData_"+varName,"");
02979 RooDataSet *projPlotData(0);
02980 if (""!=ppdStr) projPlotData=_datasets->getData(ppdStr);
02981 if(projPlotData) chkBlind(projPlotData->GetName());
02982 if(!projPlotData) projPlotData=projData;
02983 cout<<" Using dataset "<<projPlotData->GetName()<<endl;
02984 RooArgSet projVarSet;
02985
02986 RooArgSet depVarSet(_fObsSet);
02987
02988 depVarSet.remove(*theVar, kFALSE, kTRUE);
02989 depVarSet.remove(_conditionalObs, kFALSE, kTRUE);
02990
02991
02992
02993
02994
02995
02996 TString LRatioCut=readConfStr("projLRatioCut_"+varName, "-1", _runSec);
02997 if ("-1"==LRatioCut) LRatioCut=readConfStr("projLRatioCut", "-1", _runSec);
02998
02999 TString findOptCut=readConfStr("projFindOptimCut_"+varName,"notSet",_runSec);
03000 if ("notSet"==findOptCut)
03001 findOptCut=readConfStr("projFindOptimCut", "no", _runSec);
03002
03003 RooDataSet* sliceData = 0;
03004 bool delDS = false;
03005
03006 if ( LRatioCut!="-1" || findOptCut!="no" ) {
03007 const RooAbsReal *theSProj=_theSPdf->createPlotProjection(depVarSet, projVarSet);
03008 const RooAbsReal *theBProj=_theBPdf->createPlotProjection(depVarSet, projVarSet);
03009 RooArgList projPdfList(*theSProj, *theBProj, "Projected Pdf List");
03010 TString formulaStr="@0/(@0+@1)";
03011 TString funcName="lRatioFunc_"+varName;
03012 RooFormulaVar *lRatioFunc=new
03013 RooFormulaVar(funcName, "LRatio Func", formulaStr, projPdfList);
03014
03015 projPlotData=getLLRDataset(projPlotData, lRatioFunc);
03016 delDS = true;
03017
03018 if ("no"!=findOptCut) {
03019 cout<<"Finding optimal cut"<<endl;
03020
03021 TString projOptStep=readConfStr("projOptimStep_"+varName, "-1", _runSec);
03022 if ("-1"==projOptStep)
03023 projOptStep=readConfStr("projOptimStep", ".005", _runSec);
03024 Int_t sStep=(Int_t)(.5+1./atof(projOptStep));
03025 if (sStep<=0) sStep=10;
03026 cout<<"Searching step: "<<1./sStep<<", ie "<<sStep<<" steps"<<endl;
03027
03028 TString ppdCut="1";
03029 rarStrParser cutParser=ppdStr;
03030 if (cutParser.nArgs()>=2) ppdCut=cutParser[1];
03031
03032 TString optRangeCut=readConfStr("projOptimRange_"+varName, "-1", _runSec);
03033 if ("-1"==optRangeCut)
03034 optRangeCut=readConfStr("projOptimRange", "1", _runSec);
03035 if (("1"!=optRangeCut)||("1"!=ppdCut)) {
03036
03037 optRangeCut.ReplaceAll("\"", "");
03038 optRangeCut="("+optRangeCut+")&&("+ppdCut+")";
03039 cout<<"With cuts "<<optRangeCut<<endl;
03040 }
03041
03042 cout<<" Reducing "<<projPlotData->GetName()<<endl;
03043 RooDataSet *theData=(RooDataSet*)projPlotData->reduce(optRangeCut);
03044
03045 TH1F *nTotHist=new TH1F("nTotHist_"+varName,"nTotHist_"+varName,sStep,0,1);
03046 cout<<" Creating nTotHist_"+varName<<endl;
03047 for (Int_t idxEvt=0; idxEvt<theData->numEntries(); idxEvt++) {
03048 RooArgSet *dEntry=(RooArgSet*)theData->get(idxEvt);
03049 nTotHist->Fill(((RooRealVar*)dEntry->find(funcName))->getVal(),
03050 theData->weight());
03051 }
03052 delete theData;
03053 nTotHist->Print("v");
03054 plotList.Add(nTotHist);
03055
03056 rarStrParser optDataStrParser=readConfStr("projOptimData", "", _runSec);
03057
03058 Int_t projDataSize=atoi(readConfStr("projOptimDataLimit","40000",_runSec));
03059 if (projDataSize<=1000) projDataSize=1000;
03060
03061 TList hList;
03062 RooArgSet fullObsWOEVars(*_fullObs);
03063 fullObsWOEVars.remove(_protDataEVars, kFALSE, kTRUE);
03064 const char *dummy("");
03065 for (Int_t i=0; i<_sCompList.GetSize(); i++) {
03066 rarBasePdf *rarPdf=(rarBasePdf*)_sCompList.At(i);
03067 RooDataSet *theData=rarPdf->getData(dummy);
03068 if (optDataStrParser.nArgs()>0) {
03069 TString theDataName=optDataStrParser[0];
03070 optDataStrParser.Remove();
03071 theData=_datasets->getData(theDataName);
03072 if (!theData) {
03073 cout<<"Can not find dataset named "<<theDataName<<endl
03074 <<"Use the default dataset"<<endl;
03075 theData=rarPdf->getData(dummy);
03076 }
03077 }
03078 Bool_t reducedData(kFALSE);
03079
03080 if ((theData->numEntries()>projDataSize)||(_protDataEVars.getSize()>0)) {
03081 cout<<" Reducing "<<theData->GetName()<<endl;
03082 Int_t nEvt=theData->numEntries();
03083 if (nEvt>projDataSize) nEvt=projDataSize;
03084 theData=(RooDataSet*)
03085 theData->reduce(SelectVars(fullObsWOEVars),EventRange(0, nEvt-1));
03086 reducedData=kTRUE;
03087 }
03088
03089 if (_protDataEVars.getSize()>0) {
03090 RooDataSet protEData("protEData", "protEData", _protDataEVars);
03091 Int_t nEvt=theData->numEntries();
03092 Int_t nProtEvt=projPlotData->numEntries();
03093 cout<<" Getting protDataEVars from "<<projPlotData->GetName()<<endl;
03094 for (Int_t j=0; j<nEvt; j++) {
03095 Int_t I=RooRandom::randomGenerator()->Integer(nProtEvt);
03096 protEData.add(*projPlotData->get(I));
03097 }
03098 cout<<" Merging protDataEVars to "<<theData->GetName()<<endl;
03099 theData->merge(&protEData);
03100 }
03101 cout<<" Adding LRatio to "<<theData->GetName()<<endl;
03102 RooDataSet *tmpData=theData;
03103 RooDataSet *theDataNew=getLLRDataset(theData, lRatioFunc);
03104 theData=theDataNew;
03105 Double_t sigScale=theData->numEntries();
03106 theData=(RooDataSet*)theData->reduce(optRangeCut);
03107 if (reducedData) delete tmpData;
03108 sigScale=theData->numEntries()/sigScale;
03109 TString hName=Form("hist_%s_%s", rarPdf->GetName(), varName.Data());
03110 TH1F *h=new TH1F(hName, hName, sStep, 0, 1);
03111 cout<<" Creating histogram "<<h->GetName()<<endl;
03112 for (Int_t idxEvt=0; idxEvt<theData->numEntries(); idxEvt++) {
03113 RooArgSet *dEntry=(RooArgSet*)theData->get(idxEvt);
03114 h->Fill(((RooRealVar*)dEntry->find(funcName))->getVal(),
03115 theData->weight());
03116 }
03117 plotList.Add(h);
03118 h->Scale(sigScale/h->GetEntries());
03119 h->Print("v");
03120 hList.Add(h);
03121 delete theData;
03122 delete theDataNew;
03123 }
03124
03125 TString optimFormulaStr=
03126 readConfStr("projOptimFormula_"+varName, "notSet", _runSec);
03127 if ("notSet"==optimFormulaStr)
03128 optimFormulaStr=readConfStr("projOptimFormula", "notSet", _runSec);
03129 if ("notSet"==optimFormulaStr) optimFormulaStr="@0*@0/@1";
03130 optimFormulaStr.ReplaceAll("\"", "");
03131 RooRealVar NprojVar("NprojVar", "NprojVar", 0);
03132 RooRealVar NtotalVar("NtotalVar", "NtotalVar", 1);
03133 RooFormulaVar optimFormula("optimFormula", optimFormulaStr,
03134 RooArgList(NprojVar, NtotalVar));
03135 cout<<" Optimization formula: "<<optimFormulaStr<<endl;
03136 Double_t bestRatio(0.), bestCut(0.);
03137 for (Int_t i=0; i<sStep+1; i++) {
03138
03139
03140 NtotalVar.setVal(nTotHist->Integral(i, sStep));
03141 if (NtotalVar.getVal()<=0) continue;
03142 Double_t nSig(0.);
03143 for (Int_t j=0; j<_sCoefList.getSize(); j++) {
03144 nSig+=((RooAbsReal &)_sCoefList[j]).getVal()*
03145 ((TH1F*)hList.At(j))->Integral(i, sStep);
03146 }
03147 NprojVar.setVal(nSig);
03148
03149
03150 Double_t ratio=optimFormula.getVal();
03151 cout<<i<<" nSig:"<<NprojVar.getVal()<<" nTot:"<<NtotalVar.getVal()
03152 <<" ratio:"<<ratio<<endl;
03153 if (ratio>bestRatio) {
03154 bestRatio=ratio;
03155 bestCut=i;
03156 }
03157 }
03158 bestCut=bestCut/((Double_t)sStep);
03159 LRatioCut=Form("%f", bestCut);
03160 cout<<"optimal cut @ "<<bestCut<<" with n^2/ntot="<<bestRatio<<endl;
03161 }
03162 LRatioCut=funcName+">"+LRatioCut;
03163
03164 sliceData=(RooDataSet*)projPlotData->reduce(LRatioCut);
03165 } else {
03166 sliceData = (RooDataSet*)projPlotData->Clone();
03167 }
03168
03169 Int_t nBins=atoi(readConfStr(Form("plotBins_%s", theVar->GetName()),
03170 "0", _runSec));
03171 if (nBins<=0) nBins=theVar->getBins();
03172
03173 Double_t xerrorscale = atof(readConfStr("XErrorSize", "1.0", _runSec));
03174
03175 Double_t plotMin=theVar->getMin();
03176 Double_t plotMax=theVar->getMax();
03177 RooBinning *abins=
03178 getRange(theVar, "plotRange_", plotMin, plotMax, _runSec, &nBins);
03179 abins->Print();
03180
03181 TString frameName=Form("proj_%s", theVar->GetName());
03182 TString frameTitle;
03183 if ( LRatioCut != "" ) {
03184 frameTitle =
03185 Form("projection of %s with %s",theVar->GetTitle(), LRatioCut.Data());
03186 } else {
03187 frameTitle = Form("projection of %s",theVar->GetTitle());
03188 }
03189
03190 TString projAsymPlot=readConfStr("projAsymPlot_"+varName,"notSet", _runSec);
03191 if ("notSet"==projAsymPlot)
03192 projAsymPlot=readConfStr("projAsymPlot","no", _runSec);
03193
03194 TList projDS;
03195 TList projNames;
03196 projDS.Add(sliceData);
03197 projNames.Add(new TObjString(frameName));
03198
03199 TString projPlotCat=readConfStr("projPlotCat_"+varName,"notSet", _runSec);
03200 if ("notSet"==projPlotCat)
03201 projPlotCat=readConfStr("projPlotCat","no", _runSec);
03202 if ("no"!=projPlotCat) {
03203 rarStrParser projPlotCatParser=projPlotCat;
03204 while (projPlotCatParser.nArgs()>0) {
03205 TString catName=projPlotCatParser[0];
03206 projPlotCatParser.Remove();
03207 RooAbsCategory *projCat=(RooAbsCategory *)_fullObs->find(catName);
03208 if (!projCat) {
03209 cout<<" W A R N I N G !"<<endl
03210 <<" category "<<catName<<" does not exist!"<<endl;
03211 continue;
03212 }
03213 if (catName==projAsymPlot) {
03214 cout<<" "<<catName<<" is used for asymPlot"<<endl
03215 <<" Ignored"<<endl;
03216 continue;
03217 }
03218 TIterator *cIter=projCat->typeIterator();
03219 RooCatType *cType(0);
03220 while(cType=(RooCatType*)cIter->Next()) {
03221 TString catCut=catName+"=="+catName+"::"+cType->GetName();
03222 RooDataSet *catSliceData=(RooDataSet*)sliceData->reduce(catCut);
03223 projDS.Add(catSliceData);
03224 projNames.Add(new TObjString(frameName+"_"+catCut));
03225 }
03226 delete cIter;
03227 }
03228 }
03229
03230 for (Int_t pIdx=0; pIdx<projDS.GetSize(); pIdx++) {
03231 sliceData=(RooDataSet*)projDS.At(pIdx);
03232 frameName=projNames.At(pIdx)->GetName();
03233 if ("no"==projAsymPlot) {
03234 frame=getProjPlot(theVar, plotMin, plotMax, nBins, sliceData,
03235 frameName, frameTitle, plotList, RooCmdArg(), RooCmdArg(),
03236 XErrorSize(xerrorscale));
03237 } else {
03238
03239 RooCategory *cat=(RooCategory *)_fullObs->find(projAsymPlot);
03240 if (!cat) {
03241 cout<<"Can not find category: "<<projAsymPlot<<endl;
03242 exit(-1);
03243 }
03244
03245 if (2!=cat->numTypes()) {
03246 cout<<"Asym plot is for two-type category only"<<endl;
03247 exit(-1);
03248 }
03249 TIterator *catI=cat->typeIterator();
03250 RooCatType *theCatType=(RooCatType*)catI->Next();
03251 TString catCutA=
03252 projAsymPlot+"=="+projAsymPlot+"::"+theCatType->GetName();
03253 Int_t catIdxA=theCatType->getVal();
03254 theCatType=(RooCatType*)catI->Next();
03255 TString catCutB=
03256 projAsymPlot+"=="+projAsymPlot+"::"+theCatType->GetName();
03257 Int_t catIdxB=theCatType->getVal();
03258 TString catCutMinus=catIdxA<catIdxB?catCutA:catCutB;
03259 TString catCutPlus =catIdxA>catIdxB?catCutA:catCutB;
03260 RooDataSet *catSliceData=(RooDataSet*)sliceData->reduce(catCutMinus);
03261 cout<<"Dataset (-) with cut: "<<catCutMinus<<endl;
03262 catSliceData->Print();
03263 RooPlot *frameM=
03264 getProjPlot(theVar, plotMin, plotMax, nBins, catSliceData,
03265 frameName+"_Minus", frameTitle+" Minus", plotList,
03266 RooCmdArg(), RooCmdArg(), XErrorSize(xerrorscale));
03267 delete catSliceData;
03268 catSliceData=(RooDataSet*)sliceData->reduce(catCutPlus);
03269 cout<<"Dataset (+) with cut: "<<catCutPlus<<endl;
03270 catSliceData->Print();
03271 RooPlot *frameP=
03272 getProjPlot(theVar, plotMin, plotMax, nBins, catSliceData,
03273 frameName+"_Plus", frameTitle+" Plus", plotList,
03274 RooCmdArg(), RooCmdArg(), XErrorSize(xerrorscale));
03275 delete catSliceData;
03276 delete catI;
03277
03278 frame=getProjPlot(theVar, plotMin, plotMax, nBins, sliceData,
03279 frameName+"_Asym", frameTitle+" Asym", plotList,
03280 Asymmetry(*cat), Binning(*abins),
03281 XErrorSize(xerrorscale), frameM, frameP);
03282 }
03283 sliceData->Print("v");
03284 }
03285 projDS.Delete();
03286 projNames.Delete();
03287
03288
03289 TString saveDS=readConfStrCnA("projPlotSaveLLR", "no");
03290 if ("no"==saveDS) {
03291 if ( delDS ) delete projPlotData;
03292 } else {
03293 TString dsName=projPlotData->GetName();
03294 dsName+="_";
03295 dsName+=theVar->GetName();
03296 projPlotData->SetName(dsName);
03297
03298 TTree *theDSTree = createTreeFromDataset(projPlotData, kFALSE);
03299 plotList.Add(theDSTree);
03300 }
03301
03302 cout<<endl<<"lRatio Cut: "<<LRatioCut<<endl;
03303 return frame;
03304 }
03305
03325 RooPlot *rarMLFitter::getProjPlot(RooRealVar *theVar, Double_t plotMin,
03326 Double_t plotMax, Int_t nBins,
03327 RooDataSet *sliceData, TString frameName,
03328 TString frameTitle, TList &plotList,
03329 const RooCmdArg &asymCat,
03330 const RooCmdArg &nuBins,
03331 const RooCmdArg &xerrorscale,
03332 RooPlot *frameM, RooPlot *frameP)
03333 {
03334
03335 setColLimits(sliceData);
03336 RooPlot *frame(0);
03337 frame=theVar->frame(plotMin, plotMax, nBins);
03338 RooLinkedList plotOpts;
03339 plotOpts.Add((TObject*)&asymCat);
03340 plotOpts.Add((TObject*)&nuBins);
03341 plotOpts.Add((TObject*)&xerrorscale);
03342 RooCmdArg datErrArg = DataError(RooAbsData::SumW2);
03343 if (sliceData->isWeighted())
03344 plotOpts.Add((TObject*)&datErrArg);
03345 sliceData->plotOn(frame, plotOpts);
03346 Int_t lineStyle=1;
03347 Int_t lineWidth=2;
03348 Int_t lineColor=kBlue;
03349 if (!(frameM&&frameP)) {
03350
03351 _thePdf->plotOn(frame, ProjWData(*sliceData),asymCat,
03352 LineStyle(lineStyle), LineColor(lineColor));
03353 for (Int_t i=0; i<_fCompList.GetSize(); i++) {
03354 if (lineStyle>2) lineStyle=2;
03355 else lineStyle=4;
03356 rarBasePdf *compPdf=(rarBasePdf*)_fCompList.At(i);
03357 TString pdfName=compPdf->getPdf()->GetName();
03358 TString compStr=pdfName+","+pdfName+"_*";
03359 cout<<"Comps to proj: "<<compStr<<endl;
03360 _thePdf->plotOn(frame, Components(compStr), ProjWData(*sliceData),
03361 asymCat, LineWidth(lineWidth), LineStyle(lineStyle),
03362 LineColor(getColor(i)), Name(pdfName));
03363 }
03364 } else {
03365 Int_t nCurves=(Int_t)frameP->numItems();
03366 for (Int_t i=1; i<nCurves; i++) {
03367 RooCurve *asymCurve(0);
03368 RooCurve *B0Curve(0), *B0barCurve(0);
03369 B0Curve= (RooCurve*) frameP->getObject(i);
03370 B0barCurve = (RooCurve*) frameM->getObject(i);
03371 if (B0Curve&&B0barCurve) {
03372 TString curveName=Form("asymCurve_%s", B0Curve->GetName());
03373 Int_t lineWidth=B0Curve->GetLineWidth();
03374 Int_t lineStyle=B0Curve->GetLineStyle();
03375 Int_t lineColor=B0Curve->GetLineColor();
03376 asymCurve=combCurves(B0Curve, B0barCurve, "(@0-@1)/(@0+@1)", "@0+@1>0",
03377 curveName, lineWidth, lineStyle, lineColor);
03378 }
03379 if (asymCurve)
03380 frame->addPlotable(asymCurve);
03381 }
03382 }
03383
03384 frame->SetNameTitle(frameName, frameTitle);
03385
03386 frame->SetTitleSize(0.05);
03387
03388 plotList.Add(frame);
03389
03390 setColLimits(sliceData, kFALSE);
03391 return frame;
03392 }
03393
03407 RooCurve *rarMLFitter::combCurves(RooCurve *crv1, RooCurve *crv2,
03408 const char *formula, const char *valid,
03409 const char *crvName, Int_t lineWidth,
03410 Int_t lineStyle, Int_t lineColor)
03411 {
03412 RooCurve* newCurve = new RooCurve();
03413 if (0 != crvName) newCurve->SetName(crvName);
03414 newCurve->SetLineWidth(lineWidth);
03415 newCurve->SetLineStyle(lineStyle);
03416 newCurve->SetLineColor(lineColor);
03417 Double_t x1 = -1.E30;
03418 Double_t x2 = -1.E30;
03419 Double_t y1, y2;
03420 Int_t n1 = crv1->GetN();
03421 Int_t n2 = crv2->GetN();
03422 RooRealVar y1RV("y1RV", "", 0.), y2RV("y2RV", "", 0.);
03423 RooFormulaVar yComb("yComb", "", formula, RooArgSet(y1RV, y2RV));
03424 RooFormulaVar testOK("testOK", "", valid, RooArgSet(y1RV, y2RV));
03425
03426 bool needPt = true;
03427 int i(0), j(0);
03428 while (i < n1-1 && j < n2-1) {
03429 while (i < n1-1 && (needPt || x1 < x2)) {
03430 needPt = false;
03431 crv1->GetPoint(++i, x1, y1);
03432 y1RV.setVal(y1);
03433 }
03434 while (j < n2-1 && (needPt || x2 < x1)) {
03435 needPt = false;
03436 crv2->GetPoint(++j, x2, y2);
03437 y2RV.setVal(y2);
03438 }
03439 if (x1 == x2 && 0 != testOK.getVal()) {
03440
03441
03442 newCurve->addPoint(x1, yComb.getVal());
03443 }
03444 needPt = true;
03445 }
03446
03447
03448 return newCurve;
03449 }
03450
03456 void rarMLFitter::scanVarShiftToNorm(RooArgList scanVars, TArrayD &scanVarDiff)
03457 {
03458 if ("no"==readConfStr("scanVarShiftToNorm", "no", _runSec)) return;
03459 for (Int_t i=0; i<scanVars.getSize(); i++) {
03460 RooRealVar *theVar=dynamic_cast<RooRealVar*>(scanVars.at(i));
03461 if (!theVar) continue;
03462 theVar->setVal(theVar->getVal()+scanVarDiff[i]);
03463
03464 }
03465 }
03466
03476 RooPlot *rarMLFitter::doScanPlot(TList &plotList)
03477 {
03478 cout<<endl<<" In rarMLFitter doScanPlot for "<<GetName()<<endl;
03479 RooPlot *frame(0);
03480 RooCurve *curve(0);
03481
03482 TString fitOption=readConfStr("scanPlotFitOption", "qemhr", _runSec);
03483 cout<<"scanPlot fit option: \""<<fitOption<<"\""<<endl;
03484
03485
03486 Bool_t scanPlotExtended = fitOption.Contains("e");
03487 Bool_t scanPlotMinos = fitOption.Contains("m");
03488 Bool_t scanPlotHesse = fitOption.Contains("h");
03489 Bool_t scanPlotVerbose = !fitOption.Contains("q");
03490 Bool_t scanPlotSave = fitOption.Contains("r");
03491
03492
03493 RooDataSet *scanPlotData=
03494 _datasets->getData(readConfStrCnA("scanPlotData", "notSet"));
03495 if(!scanPlotData) {
03496 cout<<endl<<" Please specify in section ["<<_runSec<<"]"<<endl
03497 <<" the dataset you want to scan with config "<<endl<<endl
03498 <<" scanPlotData = <datasetName>"<<endl
03499 <<endl<<" and ru-run your job"<<endl;
03500 exit(-1);
03501 }
03502 chkBlind(scanPlotData->GetName());
03503 cout<<" Using dataset "<<scanPlotData->GetName()<<endl;
03504
03505 RooArgSet fullParams(*_thePdf->getParameters(scanPlotData));
03506
03507 string paramSStr0;
03508 writeToStr(fullParams, paramSStr0);
03509
03510 RooArgSet scanVars;
03511 rarStrParser scanVarsParser=readConfStr("scanVars", "notSet", _runSec);
03512 while (scanVarsParser.nArgs()>0) {
03513 TString varName=scanVarsParser[0];
03514 scanVarsParser.Remove();
03515 RooRealVar *theVar=(RooRealVar*)fullParams.find(varName);
03516 if (!theVar) continue;
03517 scanVars.add(*theVar);
03518
03519 Double_t min=theVar->getMin();
03520 Double_t max=theVar->getMax();
03521
03522 if (isNumber(scanVarsParser[0])) {
03523 min=atof(scanVarsParser[0]);
03524 scanVarsParser.Remove();
03525 }
03526 if (isNumber(scanVarsParser[0])) {
03527 max=atof(scanVarsParser[0]);
03528 scanVarsParser.Remove();
03529 }
03530 if (min>max) {
03531 Double_t v=min;
03532 min=max;
03533 max=v;
03534 }
03535 theVar->setRange(min,max);
03536 theVar->setConstant();
03537 }
03538 if (scanVars.getSize()<1) {
03539 cout<<" Please specify var to scan with config scanVars"<<endl;
03540 exit(-1);
03541 }
03542
03543 if (scanVars.getSize()==1) {
03544 RooRealVar *theVar=(RooRealVar*)RooArgList(scanVars).at(0);
03545 frame=theVar->frame(theVar->getMin(), theVar->getMax(), 100);
03546 frame->SetNameTitle(Form("NLLScanPlot_%s", theVar->GetName()),
03547 Form("NLLScanPlot_%s", theVar->GetName()));
03548 frame->SetYTitle("-2 ln (L/L_{0})");
03549 curve=new RooCurve();
03550 curve->SetName("NLL_curve");
03551 frame->addPlotable(curve);
03552 frame->SetMinimum(0);
03553 plotList.Add(frame);
03554 }
03555
03556 string paramSStr;
03557 writeToStr(fullParams, paramSStr);
03558
03559 Int_t nPoints=atoi(readConfStr("nScanPoints", "100", _runSec));
03560 if (nPoints<1) nPoints=1;
03561 RooRealVar NLL("NLL", "NLL", 0);
03562 RooArgSet scanSet(scanVars);
03563 scanSet.add(NLL);
03564
03565 RooDataSet *theDS=new RooDataSet("scanDS", "scanDS", scanSet);
03566
03567 TArrayD scanVarDiff(scanVars.getSize());
03568 RooArgList scanList(scanVars);
03569 for (Int_t i=0; i<scanList.getSize(); i++) {
03570 scanVarDiff[i]=((RooAbsReal&)scanList[i]).getVal();
03571 }
03572
03573 Double_t mNLL(0), maxNLL(0);
03574 scanVars.setAttribAll("Constant", kFALSE);
03575 cout<<" Refit to find mins for scanPlot"<<endl;
03576 RooFitResult *fr=_thePdf->fitTo(*scanPlotData, ConditionalObservables(_conditionalObs),
03577 Save(scanPlotSave),Extended(scanPlotExtended),
03578 Verbose(scanPlotVerbose), Hesse(scanPlotHesse),
03579 Minos(scanPlotMinos));
03580
03581
03582
03583
03584 if (fr->status()) {
03585 cout<<" Fit status for inital fit: "<<fr->status()<<endl;
03586 exit(-1);
03587 }
03588 mNLL=2*fr->minNll();
03589
03590 for (Int_t i=0; i<scanList.getSize(); i++) {
03591 scanVarDiff[i]-=((RooAbsReal&)scanList[i]).getVal();
03592 cout<<" scanVarDiff["<<scanList[i].GetName()<<"]="<<scanVarDiff[i]<<endl;
03593 }
03594
03595 Double_t theVarNormVal=((RooRealVar*)RooArgList(scanVars).at(0))->getVal();
03596 Double_t theVarNormNLL=2*fr->minNll();
03597
03598 {
03599 scanVars.setAttribAll("Constant");
03600 RooFitResult *fr=_thePdf->fitTo(*scanPlotData, ConditionalObservables(_conditionalObs),
03601 Save(scanPlotSave),Extended(scanPlotExtended),
03602 Verbose(scanPlotVerbose), Hesse(scanPlotHesse),
03603 Minos(scanPlotMinos));
03604
03605
03606
03607 if (fr->status()) {
03608 cout<<" Fit status for fit: "<<fr->status()<<endl;
03609 exit(-1);
03610 }
03611
03612 scanVarShiftToNorm(scanVars, scanVarDiff);
03613
03614 theVarNormNLL=2*fr->minNll();
03615 theVarNormVal=((RooRealVar*)RooArgList(scanVars).at(0))->getVal();
03616
03617 NLL.setVal(theVarNormNLL-mNLL);
03618 theDS->add(scanSet);
03619 }
03620
03621 Int_t nSegs=atoi(readConfStr("nScanSegments", "1", _runSec));
03622 if (nSegs<1) {
03623 cout<<" nScanSegments = "<<nSegs<<endl
03624 <<" reset to 1"<<endl;
03625 nSegs=1;
03626 }
03627 Int_t segIdx=_toyID%nSegs;
03628 for(Int_t i=0; i<nPoints; i++) {
03629
03630 readFromStr(fullParams, paramSStr);
03631
03632 RooArgList scanList(scanVars);
03633 for(Int_t j=0; j<scanList.getSize(); j++) {
03634 RooRealVar *theVar=(RooRealVar*)scanList.at(j);
03635 Double_t min=theVar->getMin();
03636 Double_t max=theVar->getMax();
03637 if (scanList.getSize()>1) {
03638 theVar->setVal(min+(max-min)*RooRandom::randomGenerator()->Uniform());
03639 } else {
03640
03641 Double_t segSize=(max-min)/nSegs;
03642
03643 min=min+segSize*segIdx;
03644 max=min+segSize;
03645 theVar->setVal(min+(max-min)*i/nPoints);
03646 }
03647 cout<<" Set scan var "<<theVar->GetName()<<" to "
03648 <<theVar->getVal()<<endl;
03649 }
03650 RooFitResult *fr=_thePdf->fitTo(*scanPlotData, ConditionalObservables(_conditionalObs),
03651 Save(scanPlotSave),Extended(scanPlotExtended),
03652 Verbose(scanPlotVerbose), Hesse(scanPlotHesse),
03653 Minos(scanPlotMinos));
03654
03655
03656
03657 if (fr->status()) {
03658 cout<<" Fit status for point #"<<i<<": "<<fr->status()<<endl;
03659 scanVars.Print("v");
03660 continue;
03661 }
03662
03663 scanVarShiftToNorm(scanVars, scanVarDiff);
03664
03665 NLL.setVal(2*fr->minNll()-mNLL);
03666 theDS->add(scanSet);
03667 if (curve) {
03668 RooRealVar *theVar=(RooRealVar*)RooArgList(scanVars).at(0);
03669 Double_t x=theVar->getVal();
03670 Double_t min=theVar->getMin();
03671 Double_t max=theVar->getMax();
03672
03673 Double_t segSize=(max-min)/nSegs;
03674
03675 min=min+segSize*segIdx;
03676 max=min+segSize;
03677 curve->addPoint(x, NLL.getVal());
03678 if ((x<theVarNormVal)&&(theVarNormVal<x+(max-min)/nPoints)) {
03679 curve->addPoint(theVarNormVal, theVarNormNLL-mNLL);
03680 }
03681 if (NLL.getVal()>maxNLL) maxNLL=NLL.getVal();
03682 }
03683 }
03684
03685 TTree *theDSTree = createTreeFromDataset(theDS, kFALSE);
03686 plotList.Add(theDSTree);
03687
03688 if (frame) frame->SetMaximum(maxNLL);
03689
03690 if (scanVars.getSize()==1) {
03691
03692 TString unCorrStr=readConfStr("scanUnCorrErr", "no", _runSec);
03693 if (!unCorrStr.BeginsWith("no")) {
03694 cout<<"scanUnCorrErr obsoleted, please use combine.C"<<endl
03695 <<"remove the config to re-run the job"<<endl;
03696 exit(-1);
03697 Double_t ucErr=atof(unCorrStr);
03698 if (ucErr<=0) {
03699 cout<<" Uncorrelated error should be > 0 ("<<unCorrStr<<")"<<endl;
03700 exit(-1);
03701 }
03702
03703 RooCurve *uccurve=frame->getCurve("NLL_curve");
03704 if (!uccurve) {
03705 cout<<" Can not find curve: NLL_curve"<<endl;
03706 exit(-1);
03707 }
03708 uccurve=(RooCurve*)uccurve->Clone();
03709 uccurve->SetName("NLL_curve_unCorr");
03710 frame->addPlotable(uccurve);
03711 addErrToCurve(uccurve, ucErr, &maxNLL);
03712 frame->SetMaximum(maxNLL);
03713 }
03714
03715
03716 TString corrStr=readConfStr("scanCorrErr", "no", _runSec);
03717 if (!corrStr.BeginsWith("no")) {
03718 cout<<"scanCorrErr obsoleted, please use combine.C"<<endl
03719 <<"remove the config to re-run the job"<<endl;
03720 exit(-1);
03721 Double_t cErr=atof(corrStr);
03722 if (cErr<=0) {
03723 cout<<" Correlated error should be > 0 ("<<corrStr<<")"<<endl;
03724 exit(-1);
03725 }
03726
03727 RooCurve *ccurve=frame->getCurve("NLL_curve_unCorr");
03728 if (!ccurve) ccurve=frame->getCurve("NLL_curve");
03729 if (!ccurve) {
03730 cout<<" Can not find curve named NLL_curve_unCorr or NLL_curve"<<endl;
03731 exit(-1);
03732 }
03733 ccurve=(RooCurve*)ccurve->Clone();
03734 ccurve->SetName("NLL_curve_Corr");
03735 frame->addPlotable(ccurve);
03736 addErrToCurve(ccurve, cErr, &maxNLL);
03737 frame->SetMaximum(maxNLL);
03738 }
03739 }
03740
03741
03742 readFromStr(fullParams, paramSStr0);
03743 return frame;
03744 }
03745
03754 void rarMLFitter::addErrToCurve(RooCurve *curve, Double_t errLo,Double_t errHi,
03755 Double_t *maxNLL)
03756 {
03757 if ((0==errLo)&&(0==errHi)) return;
03758
03759 Double_t mean(0), yMin(0), err(errLo);
03760 rarNLL nll(curve);
03761 nll.getMin(mean, yMin);
03762
03763 shiftNLLCurve(curve, 0, -yMin);
03764
03765 Int_t nPoints=curve->GetN();
03766 for(Int_t i=0; i<nPoints; i++) {
03767 Double_t x, chi2Stat;
03768 curve->GetPoint(i, x, chi2Stat);
03769 err=(x<=mean) ? errLo : errHi;
03770 if (0==err) continue;
03771 Double_t chi2Syst=(x-mean)*(x-mean)/(err*err);
03772 if (chi2Stat+chi2Syst==0.) continue;
03773 Double_t chi2Tot = chi2Stat*chi2Syst/(chi2Stat+chi2Syst);
03774 curve->SetPoint(i, x, chi2Tot);
03775 if (maxNLL) if (*maxNLL<chi2Tot) *maxNLL=chi2Tot;
03776 }
03777 }
03778
03786 void rarMLFitter::addErrToCurve(RooCurve *curve, Double_t err,Double_t *maxNLL)
03787 {
03788 addErrToCurve(curve, err, err, maxNLL);
03789 }
03790
03798 RooCurve *rarMLFitter::shiftNLLCurve(RooCurve *curve, Double_t dx, Double_t dy)
03799 {
03800 int nPoints=curve->GetN();
03801 for(Int_t i=0; i<nPoints; i++) {
03802 Double_t x(0), y(0);
03803 curve->GetPoint(i, x, y);
03804 curve->SetPoint(i, x+dx, y+dy);
03805 }
03806 return curve;
03807 }
03808
03817 RooCurve *rarMLFitter::combineNLLCurves(TList &curves, Bool_t shiftToZero,
03818 Double_t *maxNLL)
03819 {
03820
03821 set<Double_t> xs;
03822 TList nlls;
03823 for (Int_t i=0; i<curves.GetSize(); i++) {
03824 RooCurve *curve=(RooCurve*)curves.At(i);
03825 rarNLL *nll=new rarNLL(curve);
03826 nlls.Add(nll);
03827 Int_t nPoints=curve->GetN();
03828 for (Int_t j=0; j<nPoints; j++) {
03829 Double_t x, y;
03830 curve->GetPoint(j, x, y);
03831 xs.insert(x);
03832 }
03833 }
03834
03835
03836 Double_t yMax(0);
03837
03838 RooCurve *theCurve=new RooCurve();
03839 set<Double_t>::iterator the_iterator=xs.begin();
03840 while (the_iterator!=xs.end()) {
03841 Double_t y=0;
03842 for (Int_t i=0; i<curves.GetSize(); i++) {
03843 rarNLL *nll=(rarNLL*)nlls.At(i);
03844 y+=nll->getNLL(*the_iterator);
03845 }
03846 if (yMax<y) yMax=y;
03847 theCurve->addPoint(*the_iterator, y);
03848 the_iterator++;
03849 }
03850
03851 if (shiftToZero) {
03852
03853 rarNLL nll(theCurve);
03854 Double_t xMin(0), yMin(0);
03855 nll.getMin(xMin, yMin);
03856 if (yMin<-1e-05) {
03857 cout<<" Min of combined NLL is "<<yMin<<endl
03858 <<" Should be equal to or greater than 0"<<endl;
03859 }
03860
03861 shiftNLLCurve(theCurve, 0, -yMin);
03862 yMax-=yMin;
03863 }
03864
03865 if (maxNLL) if (*maxNLL<yMax) *maxNLL=yMax;
03866
03867 return theCurve;
03868 }
03869
03877 RooCurve *rarMLFitter::combineNLLCurves(TList &curves, Double_t errs[],
03878 Double_t *maxNLL)
03879 {
03880
03881 RooCurve *theCurve=combineNLLCurves(curves, kTRUE, maxNLL);
03882
03883 TList lCurves, rCurves;
03884 for(Int_t i=0; i<curves.GetSize(); i++) {
03885 RooCurve *curve=(RooCurve*)curves.At(i);
03886 RooCurve *lCurve=(RooCurve*)curve->Clone();
03887 RooCurve *rCurve=(RooCurve*)curve->Clone();
03888 shiftNLLCurve(lCurve, -fabs(errs[i]), 0);
03889 shiftNLLCurve(rCurve, +fabs(errs[i]), 0);
03890 lCurves.Add(lCurve);
03891 rCurves.Add(rCurve);
03892 }
03893 RooCurve *lCurve=combineNLLCurves(lCurves);
03894 RooCurve *rCurve=combineNLLCurves(rCurves);
03895 Double_t lx(0), rx(0), yMin(0);
03896 rarNLL lnll(lCurve), rnll(rCurve);
03897 lnll.getMin(lx, yMin);
03898 rnll.getMin(rx, yMin);
03899
03900 Double_t err=fabs(rx-lx)/2.;
03901
03902 addErrToCurve(theCurve, err, maxNLL);
03903
03904 return theCurve;
03905 }
03906
03912 Double_t rarMLFitter::getMeanErrs(RooCurve *curve,
03913 Double_t *errLo, Double_t *errHi)
03914 {
03915 Double_t theX(0), theY(0);
03916 rarNLL theNLL(curve);
03917 theNLL.getMin(theX, theY);
03918 if (fabs(theY)>1e-5) {
03919 cout<<" The NLL curve is minimized @ y="<<theY<<endl
03920 <<" (should be y=0)"<<endl;
03921 }
03922 TArrayD errs=theNLL.getX(1+theY);
03923 if (errs.GetSize()!=2) {
03924 cout<<" W A R N I N G ! !"<<endl
03925 <<" The number of points for error from NLL curve is not TWO"<<endl
03926 <<" The first and last points (if any)"
03927 <<" will be used for error calculation"<<endl
03928 <<" Your results might be problematic"<<endl
03929 <<" The points found:"<<endl;
03930 for (Int_t i=0; i<errs.GetSize(); i++) {
03931 cout<<" "<<errs[i];
03932 }
03933 cout<<endl;
03934 }
03935 if (errs.GetSize()>0) {
03936 errs[0]-=theX;
03937 errs[errs.GetSize()-1]-=theX;
03938 if (errLo) *errLo=errs[0];
03939 if (errHi) *errHi=errs[errs.GetSize()-1];
03940 }
03941 return theX;
03942 }
03943
03948 Double_t rarMLFitter::getSignf(RooCurve *curve, Double_t refVal)
03949 {
03950 rarNLL nll(curve);
03951 Double_t signf=nll.getNLL(refVal);
03952 if (signf<0) {
03953 cout<<" NLL("<<refVal<<")="<<signf<<"<0"<<endl;
03954 signf=0;
03955 }
03956 signf=sqrt(signf);
03957
03958 return signf;
03959 }
03960
03965 Double_t rarMLFitter::getUL(RooCurve *curve, Double_t CL)
03966 {
03967
03968 rarNLL theNLL(curve);
03969 Double_t theIntegral=theNLL.getLIntegral();
03970 Double_t UL=theNLL.getLIntegralInverse(theIntegral*CL);
03971 return UL;
03972 }
03973
03985 RooPlot *rarMLFitter::doContourPlot(TList &plotList)
03986 {
03987 cout<<endl<<" In rarMLFitter doContourPlot for "<<GetName()<<endl;
03988 RooPlot *frame(0);
03989
03990
03991 RooDataSet *contourPlotData=
03992 _datasets->getData(readConfStr("contourPlotData", "notSet", _runSec));
03993 if(!contourPlotData) {
03994 cout<<endl<<" Please specify in section ["<<_runSec<<"]"<<endl
03995 <<" the dataset you want to do contourPlot with config "<<endl<<endl
03996 <<" contourPlotData = <datasetName>"<<endl
03997 <<endl<<" and ru-run your job"<<endl;
03998 exit(-1);
03999 }
04000 chkBlind(contourPlotData->GetName());
04001 cout<<" Using dataset "<<contourPlotData->GetName()<<endl;
04002
04003
04004 RooArgSet fullParams(*_thePdf->getParameters(contourPlotData));
04005
04006 string paramSSaver;
04007 writeToStr(fullParams, paramSSaver);
04008
04009
04010 rarStrParser contourVarsParser=readConfStr("contourVars", "notSet", _runSec);
04011 RooRealVar *contX(0), *contY(0);
04012 Double_t xMin(0), xMax(0), yMin(0), yMax(0);
04013
04014 contX=(RooRealVar*)fullParams.find(contourVarsParser[0]);
04015 contourVarsParser.Remove();
04016 if (!contX) {
04017 cout<<"Can not find parameter "<<contourVarsParser[0]<<endl;
04018 exit(-1);
04019 }
04020 xMin=contX->getMin();
04021 xMax=contX->getMax();
04022
04023 if (isNumber(contourVarsParser[0])) {
04024 xMin=atof(contourVarsParser[0]);
04025 contourVarsParser.Remove();
04026 }
04027 if (isNumber(contourVarsParser[0])) {
04028 xMax=atof(contourVarsParser[0]);
04029 contourVarsParser.Remove();
04030 }
04031 if (xMin>xMax) {
04032 Double_t v=xMin;
04033 xMin=xMax;
04034 xMax=v;
04035 }
04036 contX->setRange(xMin,xMax);
04037
04038 contY=(RooRealVar*)fullParams.find(contourVarsParser[0]);
04039 contourVarsParser.Remove();
04040 if (!contY) {
04041 cout<<"Can not find parameter "<<contourVarsParser[0]<<endl;
04042 exit(-1);
04043 }
04044 yMin=contY->getMin();
04045 yMax=contY->getMax();
04046
04047 if (isNumber(contourVarsParser[0])) {
04048 yMin=atof(contourVarsParser[0]);
04049 contourVarsParser.Remove();
04050 }
04051 if (isNumber(contourVarsParser[0])) {
04052 yMax=atof(contourVarsParser[0]);
04053 contourVarsParser.Remove();
04054 }
04055 if (yMin>yMax) {
04056 Double_t v=yMin;
04057 yMin=yMax;
04058 yMax=v;
04059 }
04060 contY->setRange(yMin,yMax);
04061
04062 rarStrParser restParser=
04063 readConfStr("contourRestrictFloatParams", "no", _runSec);
04064 if ("no"!=restParser[0]) {
04065 Double_t dScale=2;
04066 if (isNumber(restParser[0])) {
04067 dScale=atof(restParser[0]);
04068 restParser.Remove();
04069 }
04070
04071 RooArgList fullParamList(fullParams);
04072 for (Int_t i=0; i<fullParamList.getSize(); i++) {
04073 RooRealVar *theVar=(RooRealVar*)fullParamList.at(i);
04074 if ((contX==theVar)||(contY==theVar)) continue;
04075 if (TString("RooRealVar")!=theVar->ClassName()) continue;
04076 if (theVar->isConstant()) continue;
04077 if ((!theVar->hasError()) && (!theVar->hasAsymError())) continue;
04078 Double_t min=theVar->getMin();
04079 Double_t max=theVar->getMax();
04080 Double_t myScale=dScale;
04081 Int_t myIdx=restParser.Index(theVar->GetName());
04082 if ((myIdx>=0)&&isNumber(restParser[myIdx+1]))
04083 myScale=atof(restParser[myIdx+1]);
04084 Double_t nMin(min), nMax(max);
04085 if (theVar->hasError()&&(theVar->getError()>0)) {
04086 nMin=theVar->getVal()-theVar->getError()*myScale;
04087 nMax=theVar->getVal()+theVar->getError()*myScale;
04088 }
04089 if (theVar->hasAsymError()) {
04090 if (theVar->getAsymErrorLo()<0)
04091 nMin=theVar->getVal()+theVar->getAsymErrorLo()*myScale;
04092 if (theVar->getAsymErrorHi()>0)
04093 nMax=theVar->getVal()+theVar->getAsymErrorHi()*myScale;
04094 }
04095 if (nMin>nMax) {
04096 Double_t v=nMin;
04097 nMin=nMax;
04098 nMax=v;
04099 }
04100 if (nMin>min) min=nMin;
04101 if (nMax<max) max=nMax;
04102 theVar->setRange(min, max);
04103 cout<<" Reset limits: "
04104 <<theVar->GetName()<<" ("<<min<<", "<<max<<")"<<endl;
04105 }
04106 }
04107
04108
04109 Int_t nContours=atoi(readConfStr("nContours", "2", _runSec));
04110 if (nContours<1) nContours=1;
04111 if (nContours>6) nContours=6;
04112 cout<<" Drawing "<<nContours<<" contour(s) of "
04113 <<contX->GetName()<<" vs "<<contY->GetName()<<endl<<endl;
04114 RooNLLVar nll("nll","nll",*_thePdf, *contourPlotData, kTRUE);
04115 rarMinuit min(nll);
04116 if (nContours<=1) frame = min.contour(*contX, *contY, 1, 0);
04117 else if (2==nContours) frame = min.contour(*contX, *contY, 1, 2);
04118 else if (3==nContours) frame = min.contour(*contX, *contY, 1, 2, 3);
04119 else if (4==nContours) frame = min.contour(*contX, *contY, 1, 2, 3, 4);
04120 else if (5==nContours) frame = min.contour(*contX, *contY, 1, 2, 3, 4, 5);
04121 else if (6<=nContours) frame = min.contour(*contX, *contY, 1, 2, 3, 4, 5, 6);
04122 TString myName=Form("contour_%s_%s", contX->GetName(), contY->GetName());
04123 if (frame) {
04124 frame->SetNameTitle(myName, myName);
04125 plotList.Add(frame);
04126 }
04127
04128
04129 readFromStr(fullParams, paramSSaver);
04130 return frame;
04131 }
04132
04142 RooPlot *rarMLFitter::doCombinePlot(TList &plotList)
04143 {
04144 cout<<endl<<" In rarMLFitter doCombinePlot for "<<GetName()<<endl;
04145 RooPlot *frame(0);
04146
04147
04148 Int_t ncurves = atoi(readConfStr("combineNcurves", "0", _runSec));
04149
04150 if (ncurves <= 0) {
04151 cout << "combineNcurves must be set to 1 or more curves." << endl;
04152 return (0);
04153 }
04154
04155
04156
04157 vector<Double_t> addSystErrLo, addSystErrHi, fitBias;
04158 vector<Double_t> uncorrSystErrLo, uncorrSystErrHi, corrSystErr;
04159 vector<TString> filenames, plotnames;
04160
04161
04162 Bool_t doUL = kTRUE;
04163 Double_t ConfidenceLevel(0.9);
04164 for (Int_t i=0; i < ncurves ; i++) {
04165 addSystErrLo.push_back(0.0); addSystErrHi.push_back(0.0);
04166 uncorrSystErrLo.push_back(0.0); uncorrSystErrHi.push_back(0.0);
04167 corrSystErr.push_back(0.0);
04168 fitBias.push_back(0.0);
04169 }
04170
04171 TString tempstr("");
04172
04173 tempstr = readConfStr("combineFilenames", "notSet", _runSec);
04174 if ("notSet" == tempstr) {
04175 cout << "combineFilenames must be specified." << endl;
04176 return(0);
04177 } else {
04178 rarStrParser rarFilenames = tempstr;
04179 if (rarFilenames.nArgs() < ncurves) {
04180 cout << "combineFilenames: Expected " << ncurves << " filenames, found "
04181 << rarFilenames.nArgs() << "." << endl;
04182 return(0);
04183 }
04184 for (Int_t i=0; i < rarFilenames.nArgs() ; i++) {
04185 filenames.push_back(rarFilenames[i]);
04186
04187 }
04188 }
04189
04190
04191 tempstr = readConfStr("combinePlotnames", "notSet", _runSec);
04192 if ("notSet" == tempstr) {
04193 cout << "combinePlotnames must be specified." << endl;
04194 return(0);
04195 } else {
04196 rarStrParser rarPlotnames = tempstr;
04197 if (rarPlotnames.nArgs() < ncurves) {
04198 cout << "combinePlotnames: Expected " << ncurves << " RooPlot scan plot names, found "
04199 << rarPlotnames.nArgs() << "." << endl;
04200 return(0);
04201 }
04202 for (Int_t i=0; i < rarPlotnames.nArgs() ; i++) {
04203 plotnames.push_back(rarPlotnames[i]);
04204
04205 }
04206 }
04207
04208
04209 tempstr = readConfStr("combineFitBias", "notSet", _runSec);
04210 if ("notSet" != tempstr) {
04211 rarStrParser rarFitBias = tempstr;
04212 if (rarFitBias.nArgs() < ncurves) {
04213 cout << "combineFitBias: Expected " << ncurves << " fit biases, found "
04214 << rarFitBias.nArgs() << "." << endl;
04215 return(0);
04216 }
04217 for (Int_t i=0; i < rarFitBias.nArgs() ; i++) {
04218 fitBias[i] = atof(rarFitBias[i]);
04219 }
04220 }
04221
04222
04223 tempstr = readConfStr("combineAdditive", "notSet", _runSec);
04224 if ("notSet" != tempstr) {
04225 rarStrParser rarAddErrs = tempstr;
04226 Bool_t sym = (rarAddErrs.nArgs() == ncurves);
04227 Bool_t assym = (rarAddErrs.nArgs() == (2*ncurves));
04228
04229 if (!sym && !assym) {
04230 cout << "combineAdditive: Expected either " << ncurves << " (symmetric) errors or "
04231 << ncurves*2 << " (asymmetric) errors. Found " << rarAddErrs.nArgs() << "." << endl;
04232 return (0);
04233 }
04234 for (Int_t i=0; i < ncurves ; i++) {
04235 if (assym) {
04236 addSystErrLo[i] = atof(rarAddErrs[2*i]);
04237 addSystErrHi[i] = atof(rarAddErrs[2*i+1]);
04238 } else {
04239 addSystErrHi[i] = atof(rarAddErrs[i]);
04240 addSystErrLo[i] = -1 * addSystErrHi[i];
04241 }
04242 }
04243
04244 }
04245
04246
04247 tempstr = readConfStr("combineMultiplicativeUncorrelated", "notSet", _runSec);
04248 if ("notSet" != tempstr) {
04249 rarStrParser rarSystErrs = tempstr;
04250 Bool_t sym = (rarSystErrs.nArgs() == ncurves);
04251 Bool_t assym = (rarSystErrs.nArgs() == (2*ncurves));
04252
04253 if (!sym && !assym) {
04254 cout << "combineMultiplicativeUncorrelated: Expected either "
04255 << ncurves << " (symmetric) errors or "
04256 << ncurves*2 << " (asymmetric) errors. Found "
04257 << rarSystErrs.nArgs() << "." << endl;
04258 return (0);
04259 }
04260 for (Int_t i=0; i < ncurves ; i++) {
04261
04262 if (assym) {
04263
04264 uncorrSystErrLo[i] = atof(rarSystErrs[2*i]);
04265 uncorrSystErrHi[i] = atof(rarSystErrs[2*i+1]);
04266 } else {
04267 uncorrSystErrHi[i] = atof(rarSystErrs[i]);
04268 uncorrSystErrLo[i] = -1 * uncorrSystErrHi[i];
04269 }
04270 }
04271 }
04272
04273
04274 tempstr = readConfStr("combineMultiplicativeCorrelated", "notSet", _runSec);
04275 if ("notSet" != tempstr) {
04276 rarStrParser rarSystMultErrs = tempstr;
04277 Bool_t sym = (rarSystMultErrs.nArgs() == ncurves);
04278
04279 if (!sym) {
04280 cout << "combineMultiplicativeCorrelated: Expected either " << ncurves
04281 << " errors. Found " << rarSystMultErrs.nArgs() << "." << endl;
04282 return (0);
04283 }
04284 for (Int_t i=0; i < rarSystMultErrs.nArgs() ; i++) {
04285 corrSystErr[i] = atof(rarSystMultErrs[i]);
04286 }
04287 }
04288
04289
04290 tempstr = readConfStr("combineUpperLimit", "notSet", _runSec);
04291 if ("notSet" != tempstr) {
04292 rarStrParser rarULtest = tempstr;
04293 if (rarULtest.nArgs() != 2) {
04294 cout << "combineUpperLimit: Expected <yes/no> <CL>." << endl;
04295 return (0);
04296 }
04297 if (rarULtest[0] != "yes") {doUL = kFALSE;}
04298 ConfidenceLevel = atof(rarULtest[1]);
04299 if (ConfidenceLevel <= 0.0 || ConfidenceLevel >= 100) {
04300 cout << "combineUpperLimit: Expected CL between 0 and 100%. Found "
04301 << ConfidenceLevel << "." << endl;
04302 return (0);
04303 }
04304 ConfidenceLevel *= 0.01;
04305 }
04306
04307
04308 TString xAxisTitle = "";
04309 tempstr = readConfStr("combineXaxisTitle", "notSet", _runSec);
04310 rarStrParser xtitle = tempstr;
04311 if ("notSet" != tempstr) {xAxisTitle = xtitle[0];}
04312
04313 Double_t uncorr_var(0);
04314
04315 for (Int_t j=0; j < ncurves ; j++) {
04316 cout << endl << "Curve/mode: " << j << endl;
04317 cout << "ScanPlot File: " << filenames[j] << endl;
04318 cout << "ScanPlot Name: " << plotnames[j] << endl;
04319 cout << "Additive uncorrelated Error : " << addSystErrLo[j] << "/" << addSystErrHi[j] << endl;
04320 cout << "Multiplicative correlated Error : +/- " << corrSystErr[j] << endl;
04321 cout << "Multiplicative uncorrelated Error : " << uncorrSystErrLo[j]
04322 << "/" << uncorrSystErrHi[j] << endl;
04323
04324
04325 uncorr_var = pow(uncorrSystErrHi[j],2) + pow(addSystErrHi[j],2);
04326 uncorrSystErrHi[j] = sqrt(uncorr_var);
04327
04328 uncorr_var = pow(uncorrSystErrLo[j],2) + pow(addSystErrLo[j],2);
04329 uncorrSystErrLo[j] = -1 * sqrt(uncorr_var);
04330
04331 cout << "Total Uncorrelated (add+Mult) Error : " << uncorrSystErrLo[j]
04332 << "/" << uncorrSystErrHi[j] << endl;
04333
04334 Double_t tothi = pow(uncorrSystErrHi[j],2) + pow(corrSystErr[j],2);
04335 Double_t totlo = pow(uncorrSystErrLo[j],2) + pow(corrSystErr[j],2);
04336 cout << "Total Error : " << -1 * sqrt(totlo) << "/" << sqrt(tothi) << endl;
04337 }
04338 cout << endl;
04339
04340
04341
04342 Bool_t doSignf = kTRUE;
04343 vector<Double_t> addLo, addHi, fb, unSystLo, unSystHi, corrSyst;
04344 vector<TString> fn, pn;
04345 RooPlot *thePlot(0);
04346 for (Int_t j=0; j < ncurves ; j++) {
04347
04348 fn.clear(); pn.clear(); addLo.clear(); addHi.clear();
04349 unSystLo.clear();unSystHi.clear(); corrSyst.clear(); fb.clear();
04350
04351 fn.push_back(filenames[j]); pn.push_back(plotnames[j]);
04352 addLo.push_back(addSystErrLo[j]); addHi.push_back(addSystErrHi[j]);
04353 unSystLo.push_back(uncorrSystErrLo[j]); unSystHi.push_back(uncorrSystErrHi[j]);
04354 corrSyst.push_back(corrSystErr[j]); fb.push_back(fitBias[j]);
04355
04356
04357 cout << endl << "Curve/mode: " << j << endl;
04358 thePlot = combine(1, fn, pn, fb, addLo, addHi, unSystLo, unSystHi,
04359 corrSyst, xAxisTitle, doSignf, doUL, ConfidenceLevel);
04360
04361 thePlot->SetNameTitle(Form("combinePlot_Mode%d", j),"NLL with syst+stat and stat. errs only");
04362 plotList.Add(thePlot);
04363 }
04364
04365
04366 cout << endl << "Combined curves: " << endl;
04367 thePlot = combine(ncurves, filenames, plotnames, fitBias,
04368 addSystErrLo, addSystErrHi,
04369 uncorrSystErrLo, uncorrSystErrHi, corrSystErr,
04370 xAxisTitle, doSignf, doUL, ConfidenceLevel);
04371
04372 thePlot->SetNameTitle("combinePlot","Combined curves");
04373 cout << endl;
04374 plotList.Add(thePlot);
04375
04376 return (frame);
04377 }
04378
04379
04380
04381 RooPlot *rarMLFitter::combine(Int_t nModes,
04382 const vector<TString> fileNames,
04383 const vector<TString> plotNames,
04384 const vector<Double_t> fitBias,
04385 const vector<Double_t> addSystErrLo,
04386 const vector<Double_t> addSystErrHi,
04387 const vector<Double_t> uncorrSystErrLo,
04388 const vector<Double_t> uncorrSystErrHi,
04389 const vector<Double_t> corrSystErr,
04390 const TString xAxisTitle,
04391 Bool_t doSignf, Bool_t doUL, Double_t CL)
04392 {
04393 RooPlot *thePlot(0);
04394
04395 TList curves;
04396 TList curvesWadds;
04397 TList curvesWerrs;
04398
04399
04400 TString yAxisTitle="-2 ln (L/L_{0})";
04401 Double_t xMin(0), xMax(0), yMax(0);
04402
04403
04404 for (Int_t i=0; i<nModes; i++) {
04405 TFile *f=new TFile(fileNames[i]);
04406 if (!f) {
04407 cout<<" Can not read from root file: "<<fileNames[i]<<endl;
04408 return thePlot;
04409 }
04410 RooPlot *scanPlot=(RooPlot*)f->Get(plotNames[i]);
04411 if (!scanPlot) {
04412 cout<<" Can not read in RooPlot "<<plotNames[i]
04413 <<" from "<<fileNames[i]<<endl;
04414 return thePlot;
04415 }
04416 Double_t theYmax=scanPlot->GetMaximum();
04417 if (yMax<theYmax) yMax=theYmax;
04418 TAxis *a=scanPlot->GetXaxis();
04419 Double_t theXmin=a->GetXmin();
04420 if (xMin>theXmin) xMin=theXmin;
04421 Double_t theXmax=a->GetXmax();
04422 if (xMax<theXmax) xMax=theXmax;
04423
04424
04425
04426 RooCurve *nllCurve=scanPlot->getCurve("NLL_curve");
04427 if (!nllCurve) {
04428 cout<<" Can not read NLL curve NLL_curve from RooPlot "<<plotNames[i]
04429 <<" in root file "<<fileNames[i]<<endl;
04430 return thePlot;
04431 }
04432
04433 nllCurve->SetNameTitle(Form("NLL_curve_Mode%d", i),
04434 Form("curve for sub-mode %d", i));
04435 {
04436 Double_t mean(0), yMin(0);
04437 rarNLL nll(nllCurve);
04438 nll.getMin(mean, yMin);
04439 shiftNLLCurve(nllCurve, -fitBias[i], -yMin);
04440 }
04441
04442 RooCurve *nllCurveWadd=(RooCurve*)nllCurve->Clone();
04443 RooCurve *nllCurveWerr=(RooCurve*)nllCurve->Clone();
04444
04445 addErrToCurve(nllCurveWadd, addSystErrLo[i], addSystErrHi[i]);
04446 addErrToCurve(nllCurveWerr, uncorrSystErrLo[i], uncorrSystErrHi[i]);
04447
04448 curves.Add(nllCurve);
04449 curvesWadds.Add(nllCurveWadd);
04450 curvesWerrs.Add(nllCurveWerr);
04451 }
04452
04453
04454
04455
04456
04457 RooCurve *tCurve=combineNLLCurves(curves, kTRUE, &yMax);
04458 tCurve->SetNameTitle("NLL_curve_total", "total curve w/o syst errors");
04459
04460 RooCurve *aCurve=combineNLLCurves(curvesWadds, kTRUE, &yMax);
04461 aCurve->SetNameTitle("NLL_curve_additive", "total curve w/ additve errors");
04462
04463 RooCurve *uCurve=combineNLLCurves(curvesWerrs, kTRUE, &yMax);
04464 uCurve->SetNameTitle("NLL_curve_unCorr", "total curve w/ unCorr. errors");
04465
04466
04467 Double_t corrSystErrArray[20];
04468 for (Int_t i=0; i< (Int_t) corrSystErr.size(); i++) {corrSystErrArray[i] = corrSystErr[i];}
04469 RooCurve *cCurve=combineNLLCurves(curvesWerrs, corrSystErrArray, &yMax);
04470 cCurve->SetNameTitle("NLL_curve_Corr", "total curve w/ ALL syst errors");
04471
04472 thePlot=new RooPlot(xMin, xMax, 0, yMax);
04473
04474 thePlot->SetNameTitle("combinePlot","Combined curves");
04475 thePlot->SetYTitle(yAxisTitle);
04476 thePlot->SetXTitle(xAxisTitle);
04477
04478 thePlot->addPlotable(tCurve);
04479 tCurve->SetLineWidth(2);
04480 tCurve->SetLineStyle(kDashed);
04481 thePlot->setInvisible(tCurve->GetName());
04482
04483 thePlot->addPlotable(aCurve);
04484 aCurve->SetLineWidth(2);
04485 aCurve->SetLineStyle(kDashDotted);
04486 aCurve->SetLineColor(kYellow);
04487 thePlot->setInvisible(aCurve->GetName());
04488
04489 thePlot->addPlotable(uCurve);
04490 uCurve->SetLineWidth(2);
04491 uCurve->SetLineStyle(kDashDotted);
04492 uCurve->SetLineColor(kGreen);
04493 thePlot->setInvisible(uCurve->GetName());
04494
04495 thePlot->addPlotable(cCurve);
04496 cCurve->SetLineWidth(2);
04497
04498 for (Int_t i=0; i<nModes; i++) {
04499 RooCurve *theCCurve=(RooCurve*)curves.At(i)->Clone();
04500 thePlot->addPlotable(theCCurve);
04501 theCCurve->SetLineWidth(1);
04502 theCCurve->SetLineColor(rarBasePdf::getColor(i));
04503 theCCurve->SetLineStyle(0==(i%2)?2:4);
04504 }
04505
04506 thePlot->Draw();
04507 thePlot->Print("V");
04508 cout<<endl
04509 <<" To show hidden curves (those with Option \":I\"),"
04510 <<" type, for example"<<endl
04511 <<" NLL_curve_total->SetDrawOption(\"\")"
04512 <<endl
04513 <<" To hide a curve, type, for example"<<endl
04514 <<" NLL_curve_Mode0->SetDrawOption(\"I\")"
04515 <<endl
04516 <<" Then redraw the plot: gPad->Modified()"
04517 <<endl;
04518
04519 Double_t theMean(0), systErrLo(0), systErrHi(0), statErrLo(0), statErrHi(0);
04520
04521 getMeanErrs(tCurve, &statErrLo, &statErrHi);
04522
04523 theMean=getMeanErrs(cCurve, &systErrLo, &systErrHi);
04524
04525 systErrLo=(systErrLo*systErrLo-statErrLo*statErrLo);
04526 if (systErrLo>=0) systErrLo=-sqrt(systErrLo);
04527 else {
04528 cout<<" systErrLo^2="<<systErrLo<<endl
04529 <<" Can not be true!!!! Please check!!!"<<endl;
04530 }
04531 systErrHi=(systErrHi*systErrHi-statErrHi*statErrHi);
04532 if (systErrHi>=0) systErrHi=sqrt(systErrHi);
04533 else {
04534 cout<<" systErrHi^2="<<systErrHi<<endl
04535 <<" Can not be true!!!! Please check!!!"<<endl;
04536 }
04537
04538 Double_t signfStat(0), ULStat(0), signf(0), UL(0);
04539 if (doSignf) {
04540 signfStat=getSignf(tCurve);
04541 signf=getSignf(aCurve);
04542 }
04543 if (doUL) {
04544 ULStat=getUL(tCurve, CL);
04545 UL=getUL(cCurve, CL);
04546 }
04547
04548
04549 cout<<endl;
04550 cout<<"Combined results based on "<<cCurve->GetName()
04551 <<" (and "<<tCurve->GetName()<<" for stat.):"<<endl
04552 <<" mean = "<<theMean
04553 <<" ("<<statErrLo<<", +"<<statErrHi<<")(stat)"
04554 <<" ("<<systErrLo<<", +"<<systErrHi<<")(syst)"
04555 <<endl;
04556 if (doSignf) cout<<" Signf = "<<signf<<" (wrt 0)"
04557 <<" (stat only = "<<signfStat<<")"<<endl;
04558 if (doUL) cout<<" UL = "<<UL<<" (CL="<<CL<<")"
04559 <<" (stat only = "<<ULStat<<")"<<endl;
04560
04561 return thePlot;
04562 }
04563
04574 RooPlot *rarMLFitter::doSPlot(RooRealVar *theVar, TList &plotList)
04575 {
04576 cout<<endl<<" In rarMLFitter doSPlot of "<<theVar->GetName()
04577 <<" for "<<GetName()<<endl;
04578 RooPlot *frame(0);
04579
04580 TString varName=theVar->GetName();
04581
04582 TString spdStr=readConfStrCnA("sPlotData_"+varName,"notSet");
04583 if ("notSet"==spdStr) spdStr=readConfStrCnA("sPlotData", "");
04584 RooDataSet *sPlotData(0);
04585 if (""!=spdStr) sPlotData=_datasets->getData(spdStr);
04586 if(!sPlotData) {
04587 cout<<endl<<" Please specify in section ["<<_runSec<<"]"<<endl
04588 <<" the dataset you want to project with config "<<endl<<endl
04589 <<" sPlotData = <datasetName> // and/or (see online doc)"<<endl
04590 <<" sPlotData_<obsName> = <datasetName>"<<endl
04591 <<endl<<" and ru-run your job"<<endl;
04592 exit(-1);
04593 }
04594 chkBlind(sPlotData->GetName());
04595 cout<<" Using dataset "<<sPlotData->GetName()<<endl;
04596
04597
04598 Int_t nBins=atoi(readConfStr(Form("plotBins_%s", theVar->GetName()),
04599 "0", _runSec));
04600 if (nBins<=0) nBins=theVar->getBins();
04601
04602 Double_t plotMin=theVar->getMin();
04603 Double_t plotMax=theVar->getMax();
04604 getRange(theVar, "plotRange_", plotMin, plotMax, _runSec);
04605
04606 RooArgSet ignoredObsSet(*theVar);
04607
04608 TString ignStr=readConfStr("sPlotIgnoredVars_"+varName, "notSet", _runSec);
04609 if ("notSet"==ignStr) ignStr=readConfStr("sPlotIgnoredVars", "", _runSec);
04610 rarStrParser ignoredParser=ignStr;
04611 for (Int_t i=0; i<ignoredParser.nArgs(); i++) {
04612 RooRealVar *ignored=(RooRealVar*)_fullObs->find(ignoredParser[i]);
04613 if (ignored) ignoredObsSet.add(*ignored);
04614 }
04615 cout<<" Observables not in sPlot PDF:"<<endl;
04616 ignoredObsSet.Print();
04617
04618 RooArgList thePdfsWOvar=getPdfsWOvar(ignoredObsSet);
04619
04620 if (thePdfsWOvar.getSize()<1) {
04621 cout<<" Can not build pdf for sPlot"<<endl;
04622 exit(-1);
04623 }
04624
04625 Bool_t simFit=getControlBit("SimFit");
04626
04627 if ((thePdfsWOvar.getSize()>1) && simFit) {
04628 cout<<" Currently no sPlot support for multi physCat in RooRarFit"<<endl;
04629 exit(-1);
04630 }
04631
04632 RooAbsPdf *theFullPdfWOvar=(RooAbsPdf*)thePdfsWOvar.at(0);
04633
04634 if (simFit) {
04635 RooSimPdfBuilder *simBuilder=new RooSimPdfBuilder(thePdfsWOvar);
04636 RooArgSet *simConfig=simBuilder->createProtoBuildConfig();
04637
04638 ((RooStringVar*)simConfig->find("splitCats"))->
04639 setVal(((RooStringVar*)_simConfig->find("splitCats"))->getVal());
04640 TString modelStr=((RooStringVar*)_simConfig->find("physModels"))->getVal();
04641 for (Int_t i=0; i<thePdfsWOvar.getSize(); i++) {
04642 TString pdfNameWOvar= thePdfsWOvar[i].GetName();
04643 TString pdfNameWvar=_subPdfs[i].GetName();
04644 ((RooStringVar*)simConfig->find(pdfNameWOvar))->
04645 setVal(((RooStringVar*)_simConfig->find(pdfNameWvar))->getVal());
04646 modelStr.ReplaceAll(pdfNameWvar, pdfNameWOvar);
04647 }
04648 ((RooStringVar*)simConfig->find("physModels"))->setVal(modelStr);
04649 cout<<"simPdfBuilder configs for sPlot:"<<endl;
04650 simConfig->Print("v");
04651 simBuilder->addSpecializations(_simBuilder->splitLeafList());
04652 theFullPdfWOvar=(RooAbsPdf*)
04653 simBuilder->buildPdf(*simConfig, sPlotData, 0, kTRUE);
04654
04655 }
04656 theFullPdfWOvar->Print("v");
04657
04658 RooArgSet fullParams(*theFullPdfWOvar->getParameters(sPlotData));
04659
04660 string paramSSaver;
04661 writeToStr(fullParams, paramSSaver);
04662
04663 RooArgList yields, pdfsWvar, pdfsWOvar, yield0s, pdf0sWOvar;
04664 TList rarPdfsWvar;
04665 Int_t nModel=1;
04666 if (simFit) nModel=_nComp;
04667 for (Int_t i=0; i<nModel; i++) {
04668 rarMLPdf *thisModel=(rarMLPdf*)_pdfList.At(i);
04669 TList *modelPdfList=thisModel->getPdfList();
04670 RooArgList modelCoeffList(thisModel->getSCoeffList());
04671 RooAddPdf *thisModelWOvar=(RooAddPdf*)thePdfsWOvar.at(i);
04672 RooArgList thisModelWOvarList(thisModelWOvar->pdfList());
04673 for (Int_t j=0; j<modelCoeffList.getSize(); j++) {
04674
04675 if (TString("RooRealVar")!=modelCoeffList[j].ClassName()) {
04676 cout<<" Yield "<<modelCoeffList[j].GetName()
04677 <<" must be RooRealVar"<<endl;
04678 exit(-1);
04679 }
04680 RooRealVar &theYield=(RooRealVar&)modelCoeffList[j];
04681 if (theYield.isConstant()) {
04682 cout<<" "<<theYield.GetName()<<" is constant, ignored in sPlot"<<endl;
04683 yield0s.add(theYield);
04684 } else {
04685
04686 if (theYield.getMin()>-1e6) {
04687 theYield.setMin(-1e6);
04688 cout<<" The lower limit of "<<theYield.GetName()
04689 <<" changed to -1e6"<<endl;
04690 }
04691 if (theYield.getMax()<1e6) {
04692 theYield.setMax(1e6);
04693 cout<<" The upper limit of "<<theYield.GetName()
04694 <<" changed to 1e6"<<endl;
04695 }
04696 yields.add(theYield);
04697
04698 rarBasePdf *modelCompPdf=(rarBasePdf*)modelPdfList->At(j);
04699 RooAbsPdf *thePdfWvar=modelCompPdf->getSimPdf();
04700 if (!thePdfWvar) thePdfWvar=modelCompPdf->getPdf();
04701 pdfsWvar.add(*thePdfWvar);
04702 rarPdfsWvar.Add(modelCompPdf);
04703 }
04704
04705 RooAbsPdf *thePdfWOvar=(RooAbsPdf*)thisModelWOvarList.at(j);
04706 RooAbsPdf *theSimPdfWOvar=
04707 getSimPdf((RooSimultaneous*)theFullPdfWOvar, thePdfWOvar);
04708 if (!theSimPdfWOvar) theSimPdfWOvar=thePdfWOvar;
04709 if (theYield.isConstant()) pdf0sWOvar.add(*theSimPdfWOvar);
04710 else pdfsWOvar.add(*theSimPdfWOvar);
04711 }
04712 }
04713
04714 fullParams.setAttribAll("Constant");
04715
04716 yields.setAttribAll("Constant", kFALSE);
04717
04718 TString sPlotPdfOverlay=
04719 readConfStr("sPlotPdfOverlay_"+varName,"notSet", _runSec);
04720 if ("notSet"==sPlotPdfOverlay)
04721 sPlotPdfOverlay=readConfStr("sPlotPdfOverlay", "direct", _runSec);
04722 rarStrParser pdfOverlayParser=sPlotPdfOverlay;
04723 TObjArray pdfsOverlay(yields.getSize());
04724 for (Int_t i=0; i<yields.getSize(); i++) {
04725 RooAbsPdf *OpdfI=(RooAbsPdf*)pdfsWvar.at(i);
04726 RooAbsPdf *OpdfD=((rarBasePdf*)rarPdfsWvar.At(i))->getDPdfWvar(theVar);
04727 if (!OpdfD) OpdfD=OpdfI;
04728 RooAbsPdf *theOverlayPdf=OpdfI;
04729 if ("no"==pdfOverlayParser[0]) theOverlayPdf=0;
04730 else if ("yes"!=pdfOverlayParser[0]) theOverlayPdf=OpdfD;
04731
04732 Int_t idx=pdfOverlayParser.Index(yields[i].GetName())+1;
04733 if (idx>0) {
04734 if ("no"==pdfOverlayParser[idx]) theOverlayPdf=0;
04735 else if("yes"==pdfOverlayParser[idx]) theOverlayPdf=OpdfI;
04736 else if("direct"==pdfOverlayParser[idx]) theOverlayPdf=OpdfD;
04737 else {
04738 rarBasePdf *thePdf=(rarBasePdf*)
04739 (_rarPdfs.FindObject(pdfOverlayParser[idx]));
04740 if (!thePdf)
04741 cout<<" Can not find pdf named "<<pdfOverlayParser[idx]
04742 <<" for yield "<<yields[i].GetName()
04743 <<" in config sPlotPdfOverlay"<<endl;
04744 else {
04745 RooAbsPdf *thisPdf=thePdf->getSimPdf();
04746 if (!thisPdf) thisPdf=thePdf->getPdf();
04747 theOverlayPdf=thisPdf;
04748 }
04749 }
04750 }
04751 pdfsOverlay[i]=theOverlayPdf;
04752 }
04753
04754 yields.Print();
04755 pdfsWvar.Print();
04756 pdfsWOvar.Print();
04759
04760 RooArgSet projDeps(getArgSet("sPlotNormIgnoredObs", kTRUE));
04761 projDeps.add
04762 (getArgSet(readConfStr("sPlotNormIgnoredObs","",_runSec),kFALSE));
04763
04764 projDeps.add(_conditionalObs);
04765
04766
04767
04768
04769 TString fitOption("qemhr");
04770
04771
04772 Bool_t sPlotExtended = fitOption.Contains("e");
04773 Bool_t sPlotMinos = fitOption.Contains("m");
04774 Bool_t sPlotHesse = fitOption.Contains("h");
04775 Bool_t sPlotVerbose = !fitOption.Contains("q");
04776 Bool_t sPlotSave = fitOption.Contains("r");
04777 RooFitResult *fitStat=
04778 theFullPdfWOvar->fitTo(*sPlotData, ConditionalObservables(_conditionalObs),
04779 Save(sPlotSave), Extended(sPlotExtended),
04780 Verbose(sPlotVerbose), Hesse(sPlotHesse),
04781 Minos(sPlotMinos));
04782
04783
04784
04785
04786 TString sPlotName=Form("sPlot_%s", theVar->GetName());
04787 rarSPlot mySPlot(sPlotName, sPlotName, theVar, *sPlotData, fitStat,
04788 pdfsWOvar, yields, pdf0sWOvar, yield0s, projDeps);
04789
04790 rarStrParser sPlotComps=readConfStr("sPlotComps", "all", _runSec);
04791
04792 TString sPlotAsym=readConfStr("sPlotAsym_"+varName,"notSet", _runSec);
04793 if ("notSet"==sPlotAsym)
04794 sPlotAsym=readConfStr("sPlotAsym","no", _runSec);
04795 RooCategory *sPlotAsymCat(0);
04796 if ("no"!=sPlotAsym) {
04797 sPlotAsymCat=(RooCategory *)_fullObs->find(sPlotAsym);
04798 if (!sPlotAsymCat) {
04799 cout<<"Can not find category: "<<sPlotAsym<<endl;
04800 exit(-1);
04801 }
04802
04803 if (2!=sPlotAsymCat->numTypes()) {
04804 cout<<"Asym plot is for two-type category only"<<endl;
04805 exit(-1);
04806 }
04807 }
04808
04809 TList sPlotCatList;
04810 TString sPlotCatName=readConfStr("sPlotCat_"+varName,"notSet", _runSec);
04811 if ("notSet"==sPlotCatName)
04812 sPlotCatName=readConfStr("sPlotCat","no", _runSec);
04813 if ("no"!=sPlotCatName) {
04814 rarStrParser sPlotCatParser=sPlotCatName;
04815 while (sPlotCatParser.nArgs()>0) {
04816 TString catName=sPlotCatParser[0];
04817 sPlotCatParser.Remove();
04818 RooAbsCategory *sPlotCat=(RooAbsCategory *)_fullObs->find(catName);
04819 if (!sPlotCat) {
04820 cout<<" W A R N I N G !"<<endl
04821 <<" category "<<catName<<" does not exist!"<<endl;
04822 continue;
04823 }
04824 if (catName==sPlotAsym) {
04825 cout<<" "<<catName<<" is used for sPlotAsym"<<endl
04826 <<" Ignored"<<endl;
04827 continue;
04828 }
04829 sPlotCatList.Add(sPlotCat);
04830 }
04831 }
04832 RooLinkedList plotOpts;
04833 Double_t xerrorscale = atof(readConfStr("XErrorSize", "1.0", _runSec));
04834 RooCmdArg xerrArg = XErrorSize(xerrorscale);
04835 plotOpts.Add((TObject*)&xerrArg);
04836 RooCmdArg datErrArg = DataError(RooAbsData::SumW2);
04837 RooCmdArg asymArg;
04838 if(sPlotAsymCat) {
04839 asymArg = Asymmetry(*sPlotAsymCat);
04840 plotOpts.Add((TObject*)&asymArg);
04841 }
04842 else
04843 plotOpts.Add((TObject*)&datErrArg);
04844 for (Int_t i=0; i<yields.getSize(); i++) {
04845 RooRealVar &yield=(RooRealVar&)yields[i];
04846 TString yieldName=yield.GetName();
04847 TString pdfName=pdfsWvar[i].GetName();
04848 pdfName.Remove(0,4);
04849 if (("all"!=sPlotComps[0])&&(!sPlotComps.Have(yieldName))
04850 &&(!sPlotComps.Have(pdfName))) continue;
04851 cout<<" For "<<yieldName<<" "<<varName<<endl;
04852 TString mySPlotName=sPlotName+"_"+yieldName;
04853 RooDataSet *fillData=mySPlot.fill(yield, nBins, plotMin, plotMax);
04854
04855 TList sPlotDS;
04856 TList sPlotSlice;
04857 TList sPlotNames;
04858 sPlotDS.Add(new RooDataSet(*fillData));
04859 sPlotSlice.Add(new RooDataSet(*sPlotData));
04860 sPlotNames.Add(new TObjString(mySPlotName));
04861 for( int sPlotCatIdx=0; sPlotCatIdx < sPlotCatList.GetSize();
04862 sPlotCatIdx++) {
04863 RooAbsCategory *sPlotCat=(RooAbsCategory*)sPlotCatList.At(sPlotCatIdx);
04864 TIterator *cIter=sPlotCat->typeIterator();
04865 RooCatType *cType(0);
04866 while(cType=(RooCatType*)cIter->Next()) {
04867 TString catName=sPlotCat->GetName();
04868 TString typeName=cType->GetName();
04869 TString catCut=catName+"=="+catName+"::"+typeName;
04870 RooDataSet *catFillData=(RooDataSet*)fillData->reduce(catCut);
04871 RooDataSet *catSliceData=(RooDataSet*)sPlotData->reduce(catCut);
04872 sPlotDS.Add(catFillData);
04873 sPlotSlice.Add(catSliceData);
04874 sPlotNames.Add(new TObjString(mySPlotName+"_"+typeName));
04875 }
04876 delete cIter;
04877 }
04878 for (Int_t sPlotIdx=0; sPlotIdx<sPlotDS.GetSize(); sPlotIdx++) {
04879 frame=theVar->frame(plotMin, plotMax, nBins);
04880 TString frameName=sPlotNames.At(sPlotIdx)->GetName();
04881 RooDataSet* localData = (RooDataSet*)sPlotDS.At(sPlotIdx);
04882 RooDataSet* localSlice = (RooDataSet*)sPlotSlice.At(sPlotIdx);
04883 frame->SetNameTitle(frameName,frameName);
04884
04885 if(sPlotIdx==0 &&
04886 !(localData&&"yes"!=readConfStr("sPlotHist", "no", _runSec)))
04887 frame->addTH1(mySPlot.getSPlotHist(), "E");
04888
04889 if (localData&&"yes"!=readConfStr("sPlotHist", "no", _runSec)) {
04890 if (!sPlotAsymCat)
04891 localData->plotOn(frame, plotOpts);
04892 else
04893
04894 localData->plotOn(frame, plotOpts);
04895 }
04896 RooAbsPdf *theOverlayPdf=(RooAbsPdf*)pdfsOverlay[i];
04897 if (theOverlayPdf) {
04898 if (theOverlayPdf->dependsOn(*theVar)) {
04899 cout<<" Overlay "<<theOverlayPdf->GetName()<<" on sPlot"<<endl;
04900 if (!sPlotAsymCat)
04901 theOverlayPdf->plotOn(frame, ProjWData(projDeps, *localSlice));
04902 else {
04903 TIterator *catI=sPlotAsymCat->typeIterator();
04904 RooCatType *theCatType=(RooCatType*)catI->Next();
04905 TString catCutA=
04906 sPlotAsym+"=="+sPlotAsym+"::"+theCatType->GetName();
04907 Int_t catIdxA=theCatType->getVal();
04908 theCatType=(RooCatType*)catI->Next();
04909 TString catCutB=
04910 sPlotAsym+"=="+sPlotAsym+"::"+theCatType->GetName();
04911 Int_t catIdxB=theCatType->getVal();
04912 TString catCutMinus=catIdxA<catIdxB?catCutA:catCutB;
04913 TString catCutPlus =catIdxA>catIdxB?catCutA:catCutB;
04914 RooDataSet *catLocalData=(RooDataSet*)localData->reduce(catCutMinus);
04915 RooDataSet *catLocalSlice=(RooDataSet*)localSlice->reduce(catCutMinus);
04916 cout<<"Dataset (-) with cut: "<<catCutMinus<<endl;
04917 catLocalData->Print();
04918 RooPlot *frameM=theVar->frame(plotMin, plotMax, nBins);
04919 catLocalData->plotOn(frameM, DataError(RooAbsData::SumW2));
04920 theOverlayPdf->plotOn(frameM, ProjWData(projDeps,*catLocalSlice));
04921 RooCurve *B0barCurve = (RooCurve*) frameM->getObject(1);
04922 delete catLocalData; delete catLocalSlice;
04923 catLocalData=(RooDataSet*)localData->reduce(catCutPlus);
04924 catLocalSlice=(RooDataSet*)localSlice->reduce(catCutPlus);
04925 cout<<"Dataset (+) with cut: "<<catCutPlus<<endl;
04926 catLocalData->Print();
04927 RooPlot *frameP=theVar->frame(plotMin, plotMax, nBins);
04928 catLocalData->plotOn(frameP, DataError(RooAbsData::SumW2));
04929 theOverlayPdf->plotOn(frameP, ProjWData(projDeps,*catLocalSlice));
04930 RooCurve *B0Curve = (RooCurve*) frameP->getObject(1);
04931 delete catLocalData;
04932 delete catI;
04933 TString curveName=Form("asymCurve_%s", B0Curve->GetName());
04934 RooCurve* asymCurve=combCurves(B0Curve,B0barCurve,
04935 "(@0-@1)/(@0+@1)", "@0+@1>0",
04936 curveName, 3, kSolid, kBlue);
04937 if(asymCurve)
04938 frame->addPlotable(asymCurve);
04939 }
04940 cout<<endl;
04941 } else {
04942 cout<<" Overlaying pdf "<<theOverlayPdf->GetName()
04943 <<" does not depend on "<<theVar->GetName()<<endl
04944 <<" Overlaying ignored!!"<<endl;
04945 }
04946 }
04947
04948 plotList.Add(frame);
04949 }
04950 sPlotDS.Delete();
04951 sPlotSlice.Delete();
04952 sPlotNames.Delete();
04953 if (("yes"==readConfStr("sPlotSaveSWeight", "no", _runSec))&&fillData) {
04954 TString dName=frame->GetName();
04955 dName.ReplaceAll("sPlot_", "sWeight_");
04956
04957 TTree *dataTree = createTreeFromDataset(fillData, kFALSE);
04958 plotList.Add(dataTree);
04959 } else {
04960 delete fillData;
04961 }
04962 }
04963
04964
04965 readFromStr(fullParams, paramSSaver);
04966 return frame;
04967 }
04968
04983 Bool_t rarMLFitter::initParams(TString act, RooArgSet fullParams,
04984 RooArgSet fullParamsWOI,
04985 TString readParams, TString paramFileID,
04986 TString readSecParams, Bool_t useFloatFix,
04987 RooArgSet postPdfFloatSet,
04988 RooArgSet preMLFixSet, RooArgSet preMLFloatSet)
04989 {
04990 Bool_t retVal(kTRUE);
04991
04992 TString readParamsConfig=readConfStrCnA("pre"+act+"ReadParams", readParams);
04993 if ("no"!=readParamsConfig) {
04994 TString paramFile=getParamFileName(paramFileID, readParamsConfig);
04995 cout<<endl<<"Reading in params before "<<act<<" from"<<endl
04996 <<paramFile<<endl;
04997 paramFileI(fullParamsWOI, paramFile);
04998 }
04999
05000 TString readSecParamsConfig=
05001 readConfStr("pre"+act+"ReadSecParams", readSecParams, _runSec);
05002 if ("no"!=readSecParamsConfig) {
05003 if ("yes"!=readSecParamsConfig) {
05004 cout<<endl<<"Reading in params for "<<act<<" from section \""
05005 <<readSecParamsConfig<<"\""<<endl;
05006 fullParams.readFromFile(_configFile, 0, readSecParamsConfig);
05007 }
05008 cout<<endl<<"Reading in params for "<<act<<" from section \""
05009 <<_runSec<<"\""<<endl;
05010 fullParams.readFromFile(_configFile, 0, _runSec);
05011 }
05012 if (useFloatFix) {
05013
05014 postPdfFloatSet.setAttribAll("Constant", kFALSE);
05015
05016 preMLFixSet.setAttribAll("Constant");
05017
05018 preMLFloatSet.setAttribAll("Constant", kFALSE);
05019 }
05020 cout<<endl<<" The init values of variables for "<<act<<":"<<endl;
05021 fullParams.Print("v");
05022
05023 return retVal;
05024 }
05025
05030
05031 void rarMLFitter::saveAsRootFile(const RooDataSet *ds,
05032 const TString rootfilename,
05033 Bool_t withErrors = kTRUE) {
05034
05035 cout << "Saving tree for " << ds->GetName() << endl;
05036 TFile *f = new TFile(rootfilename,"RECREATE");
05037 if (f != 0) {
05038 TTree *tree = createTreeFromDataset(ds, withErrors);
05039 if (tree != 0) {
05040 tree->Write();
05041 } else {
05042 cout << "Warning: empty tree for dataset " << ds->GetName() << endl;
05043 }
05044 f->Close();
05045
05046 }
05047 return;
05048 }
05049
05053
05054 TTree *rarMLFitter::createTreeFromDataset(const RooDataSet *ds,
05055 Bool_t withErrors = kTRUE) {
05056
05057 const Int_t NVARS = 5000;
05058 Double_t array[NVARS];
05059
05060 const Bool_t debug = kFALSE;
05061
05062
05063 const Int_t nrows = ds->numEntries();
05064
05065 if (debug) {
05066 cout << "Number of nrows " << nrows << " for ds " << ds->GetName() << endl;
05067 }
05068
05069 if (nrows < 1) {
05070 cout << "No entries in dataset " << ds->GetName() << endl;
05071 return(0);
05072 }
05073
05074
05075 const RooArgSet *names = ds->get(0);
05076 const Int_t nvars = names->getSize();
05077
05078 if (nvars > NVARS) {
05079 cout << "Too many observables )"<<nvars<<") to save. "
05080 << "Maximum is " << NVARS << endl;
05081 }
05082
05083
05084 TTree *tree = new TTree(ds->GetName(), ds->GetTitle());
05085
05086
05087 Int_t icol(0);
05088 TIterator *iter = names->createIterator();
05089
05090
05091 typedef map<TString, Int_t> MapType;
05092 MapType my_map;
05093
05094 RooRealVar *rvar(0);
05095 while ((rvar = (RooRealVar*) iter->Next())) {
05096 if (rvar->ClassName() == TString("RooRealVar")) {
05097 TString varname = rvar->GetName();
05098 my_map[varname] = icol; icol++;
05099 if (withErrors) {
05100
05101
05102
05103 if (rvar->hasError()) {
05104 if (rvar->getError() != 0) {my_map[varname+"_err"] = icol; icol++;}
05105 }
05106 if (rvar->hasAsymError()) {
05107 if (rvar->getAsymErrorLo() != 0 || rvar->getAsymErrorHi() != 0) {
05108 my_map[varname+"_aerr_lo"] = icol; icol++;
05109 my_map[varname+"_aerr_hi"] = icol; icol++;
05110 }
05111 }
05112 }
05113 if (debug) {cout << "VARR " << icol << " " << varname << endl;}
05114 }
05115 }
05116
05117
05118 RooAbsArg *var(0);
05119 iter->Reset();
05120 while ((var = (RooAbsArg*) iter->Next())) {
05121 if (var->ClassName() != TString("RooRealVar")) {
05122 TString varname = var->GetName();
05123 if (debug) {cout << "VAR C " << icol << " " << varname << endl;}
05124 my_map[varname] = icol;
05125 icol++;
05126 }
05127 }
05128
05129
05130 MapType::const_iterator it;
05131 for (it = my_map.begin(); it != my_map.end(); ++it) {
05132 TString varname = it->first;
05133 Int_t icol = it->second;
05134 tree->Branch(varname.Data(), &array[icol]);
05135 }
05136
05137
05138 for (Int_t irow=0; irow < nrows; irow++) {
05139
05140 const RooArgSet *row = ds->get(irow);
05141 icol = -1;
05142 iter->Reset() ;
05143 const Int_t idefval = -992746521;
05144 const Double_t ddefval = -992746521;
05145
05146 while ((var = (RooAbsArg*) iter->Next())) {
05147 icol++;
05148 TString varname = var->GetName();
05149 MapType::const_iterator iter1 = my_map.find(varname);
05150 if (iter1 != my_map.end()) {icol = iter1->second;}
05151
05152
05153 Int_t itemp = row->getCatIndex(varname, idefval);
05154 if (itemp != idefval) {
05155 array[icol] = (Double_t) itemp;
05156 continue;
05157 }
05158
05159 Double_t temp = row->getRealValue(varname, ddefval);
05160 if (temp != ddefval) {
05161
05162 array[icol] = temp;
05163
05164
05165 RooRealVar *rvar = dynamic_cast<RooRealVar*>(var);
05166 Double_t err(0);
05167 iter1 = my_map.find(varname+"_err");
05168 if (iter1 != my_map.end()) {
05169 icol = iter1->second;
05170 array[icol] = rvar->getError();
05171 }
05172
05173
05174 iter1 = my_map.find(varname+"_aerr_lo");
05175 if (iter1 != my_map.end()) {
05176 icol = iter1->second;
05177 if (rvar->hasAsymError()) {
05178 err = rvar->getAsymErrorLo();
05179 } else {
05180 err = 1;
05181 }
05182 array[icol] = err;
05183 }
05184
05185
05186 iter1 = my_map.find(varname+"_aerr_hi");
05187 if (iter1 != my_map.end()) {
05188 icol = iter1->second;
05189 array[icol] = rvar->getAsymErrorHi();
05190 if (rvar->hasAsymError()) {
05191 err = rvar->getAsymErrorHi();
05192 } else {
05193 err = -1;
05194 }
05195 array[icol] = err;
05196 }
05197 continue;
05198 }
05199 cout << "Unexpected type " << var->ClassName() << endl;
05200
05201 }
05202 tree->Fill();
05203 }
05204
05205 return (tree);
05206 }
05207
05227 void rarMLFitter::run()
05228 {
05229
05230 cout<<endl<<" Pre-Actions for each RooRarFitPdf"<<endl<<endl;
05231 preAction();
05232
05233 rarStrParser compCorrParser=
05234 readConfStr("computeCorrelations", "yes", getDatasets()->getVarSec());
05235
05236 for (Int_t i=0; i<getDatasets()->getDatasetList()->GetSize(); i++) {
05237 RooDataSet *theData=(RooDataSet*)getDatasets()->getDatasetList()->At(i);
05238 if (("yes"==compCorrParser[0])||compCorrParser.Have(theData->GetName()))
05239 computeCorrelations(_fObsSet, theData);
05240 }
05241
05242 cout<<endl<<endl<<" Run rarMLFitter with configs from section \""
05243 <<_runSec<<"\""<<endl<<endl;
05244
05245
05246 RooArgSet fullParams(getParams());
05247
05248 fullParams.add(*_thePdf->getParameters(_protDataset));
05249
05250 TString paramOrder=readConfStr("outParamOrder","ascend",getMasterSec());
05251 if (("ascend"==paramOrder) || ("descend"==paramOrder)) {
05252 Bool_t reverse=("descend"==paramOrder);
05253 RooArgList paramList(fullParams);
05254 fullParams.removeAll();
05255 paramList.sort(reverse);
05256 fullParams.add(paramList);
05257 }
05258
05259 RooArgSet fullParamsWOI(fullParams);
05260 fullParamsWOI.remove(getArgSet("Ignored", kTRUE, &fullParams));
05261
05262 RooArgSet postPdfFloatSet(getArgSet("postPdfFloat", kTRUE, &fullParams));
05263 postPdfFloatSet.add(getArgSet(readConfStr("postPdfFloat","", _runSec),kFALSE,
05264 &fullParams));
05265
05266 RooArgSet preMLFixSet(getArgSet("preMLFix", kTRUE, &fullParams));
05267 preMLFixSet.add(getArgSet(readConfStr("preMLFix","", _runSec),kFALSE,
05268 &fullParams));
05269
05270 RooArgSet preMLFloatSet(getArgSet("preMLFloat", kTRUE, &fullParams));
05271 preMLFloatSet.add(getArgSet(readConfStr("preMLFloat","", _runSec),kFALSE,
05272 &fullParams));
05273
05274
05275 if ("no"!=readConfStr("pdfFit", "no", _runSec)) {
05276 cout<<endl<<" Do pdfFit Action"<<endl<<endl;
05277
05278 initParams("Pdf", fullParams, fullParamsWOI, "no", "pdfFit", "no");
05279
05280 TString pdfToFit=readConfStr("pdfToFit", "", _runSec);
05281 if (""!=pdfToFit) cout<<" Pdf(s) to fit: "<<pdfToFit<<endl;
05282
05283 doPdfFit(pdfToFit);
05284
05285 TString postPdfMakePlot=readConfStr("postPdfMakePlot", "no", _runSec);
05286 if ("no"!=postPdfMakePlot) {
05287
05288 TList *plotList=new TList();
05289 doPdfPlot(*plotList, pdfToFit);
05290 TString postPdfPlotFile=getRootFileName("pdfPlot", postPdfMakePlot);
05291 TFile plotFile(postPdfPlotFile, "recreate");
05292 plotList->Print();
05293 plotList->Write();
05294 plotFile.Close();
05295 cout<<endl<<"Writing out Pdf plots after pdf fit to "
05296 <<postPdfPlotFile<<endl;
05297 }
05298
05299 TString postPdfReadSecParams="no";
05300 postPdfReadSecParams=
05301 readConfStr("postPdfReadParams", postPdfReadSecParams, _runSec);
05302 postPdfReadSecParams=
05303 readConfStr("postPdfReadSecParams", postPdfReadSecParams, _runSec);
05304 if ("no"!=postPdfReadSecParams) {
05305 if ("yes"!=postPdfReadSecParams) {
05306 cout<<endl<<"Reading in params after pdfFit from section \""
05307 <<postPdfReadSecParams<<"\""<<endl;
05308 fullParams.readFromFile(_configFile, 0, postPdfReadSecParams);
05309 }
05310 cout<<endl<<"Reading in params after pdfFit from section \""
05311 <<_runSec<<"\""<<endl;
05312 fullParams.readFromFile(_configFile, 0, _runSec);
05313 }
05314
05315
05316 cout<<endl<<"Fix params after pdfFit"<<endl;
05317 fullParams.setAttribAll("Constant");
05318
05319 TString postPdfWriteParams=readConfStrCnA("postPdfWriteParams", "no");
05320 if ("no"!=postPdfWriteParams) {
05321
05322 if (""!=pdfToFit) {
05323 fullParams.writeToStream(cout, kFALSE);
05324 cout<<endl<<"W A R N I N G !"<<endl
05325 <<"You do not fit all PDFs!"<<endl
05326 <<"No output for the params!"<<endl;
05327 } else {
05328
05329 TString postPdfParamFile=
05330 getParamFileName("pdfFit", postPdfWriteParams);
05331 cout<<endl<<"Writing out params after pdfFit to "
05332 <<postPdfParamFile<<endl;
05333 paramFileO(fullParams, postPdfParamFile);
05334 }
05335 }
05336 }
05337
05338
05339 if ("no"!=readConfStr("toyStudy", "no", _runSec)) {
05340 cout<<endl<<" Do toyStudy Action"<<endl<<endl;
05341
05342 initParams("Toy", fullParams, fullParamsWOI, "yes", "pdfFit", "yes",
05343 kTRUE, postPdfFloatSet, preMLFixSet, preMLFloatSet);
05344
05345 RooDataSet *theFitParDataSet=doToyStudy(fullParams);
05346 TString postToyWriteParams=readConfStr("postToyWriteParams","yes",_runSec);
05347 if (("no"!=postToyWriteParams) && theFitParDataSet) {
05348 TString toyRootFile=getRootFileName("toyPlot", postToyWriteParams);
05349 cout<<endl<<"Writing out param pulls, etc after toy study to "
05350 <<toyRootFile<<endl;
05351 saveAsRootFile(theFitParDataSet, toyRootFile, kTRUE);
05352 }
05353 }
05354
05355
05356 if ("no"!=readConfStr("mlFit", "no", _runSec)) {
05357 cout<<endl<<" Do mlFit Action"<<endl<<endl;
05358
05359 TString postMLSysParams=readConfStrCnA("postMLSysParams", "no");
05360 TString postMLSysVars=readConfStrCnA("postMLSysVars", "no");
05361 if((!postMLSysParams.BeginsWith("no"))&&(!postMLSysVars.BeginsWith("no")))
05362 {
05363 TString vParams="";
05364 rarStrParser paramsStrParser=postMLSysParams;
05365 while (paramsStrParser.nArgs()>0) {
05366 TString theParamName=paramsStrParser[0];
05367 paramsStrParser.Remove();
05368 TIterator* iter = fullParams.createIterator();
05369 RooRealVar *theParam(0);
05370 while(theParam=(RooRealVar*)iter->Next()) {
05371 TString theName=theParam->GetName();
05372 if ((theName==theParamName)||(theName.BeginsWith(theParamName+"_"))){
05373 vParams=vParams+" \""+theParamName+"\"";
05374 break;
05375 }
05376 }
05377 delete iter;
05378 }
05379 rarStrParser vParamsParser=vParams;
05380 for (Int_t i=0; i<vParamsParser.nArgs(); i++) {
05381 for (Int_t j=0; j<i; j++) {
05382 saveCorrCoeff(getCorrCoefName(vParamsParser[i],vParamsParser[j]),0,
05383 kTRUE);
05384 }
05385 }
05386
05387 }
05388
05389 initParams("ML", fullParams, fullParamsWOI, "yes", "pdfFit", "yes",
05390 kTRUE, postPdfFloatSet, preMLFixSet, preMLFloatSet);
05391
05392 TString mlFitOption=readConfStr("mlFitOption", "hq", _runSec);
05393
05394 mlFitOption+="er";
05395 cout<<"mlFit option: \""<<mlFitOption<<"\""<<endl;
05396
05397 stringstream o;
05398
05399 RooDataSet *mlFitData=
05400 _datasets->getData(readConfStrCnA("mlFitData", "notSet"));
05401 if(!mlFitData) {
05402 cout<<endl<<" Please specify in section ["<<_runSec<<"]"<<endl
05403 <<" the dataset you want to fit on with config "<<endl<<endl
05404 <<" mlFitData = <datasetName>"<<endl<<endl
05405 <<" and ru-run your job"<<endl;
05406 exit(-1);
05407 }
05408
05409 chkBlind(mlFitData->GetName());
05410 cout<<" Using dataset "<<mlFitData->GetName()<<" in mlAction"<<endl;
05411
05412 Int_t mlFitNumCPU=atoi(readConfStr("useNumCPU", "1", getMasterSec()));
05413
05414
05415 TList *plotList=new TList();
05416
05417 RooFitResult *fitResult=doMLFit(mlFitData, mlFitOption, o, mlFitNumCPU);
05418
05419
05420 TString postMLSignf=readConfStrCnA("postMLSignf", "no");
05421 if (!postMLSignf.BeginsWith("no"))
05422 doSignf(mlFitData, postMLSignf, fitResult, fullParams, o);
05423
05424 if ((!postMLSysParams.BeginsWith("no"))&&(!postMLSysVars.BeginsWith("no")))
05425 doSysStudy(mlFitData, postMLSysParams, postMLSysVars, fullParams, o);
05426
05427 TString postMLGOFChisq=readConfStrCnA("postMLGOFChisq", "no");
05428 if (!postMLGOFChisq.BeginsWith("no"))
05429 doGOFChisq(mlFitData, o, plotList);
05430
05431 TString mlPlotFile=readConfStr("mlPlotFile", "default", _runSec);
05432 mlPlotFile=getRootFileName("mlPlot", mlPlotFile);
05433 TFile plotFile(mlPlotFile,"recreate");
05434 plotList->Print();
05435 plotList->Write();
05436 TString postMLWriteFitResult=readConfStrCnA("postMLWriteFitResult", "no");
05437 if("no"!=postMLWriteFitResult) {
05438
05439 fitResult->Write();
05440 }
05441 plotFile.Close();
05442 cout<<"Writing out ml plots to "<<mlPlotFile<<endl;
05443
05444 cout<<o.str()<<endl;
05445
05446 TString postMLWriteParams=readConfStrCnA("postMLWriteParams", "yes");
05447 if ("no"!=postMLWriteParams) {
05448 TString postMLParamFile=getParamFileName("mlFit", postMLWriteParams);
05449 cout<<endl<<"Writing out params after mlFit to "<<postMLParamFile<<endl;
05450
05451 fullParams.setAttribAll("Constant");
05452 paramFileO(fullParams, postMLParamFile);
05453 }
05454 }
05455
05456
05457 if ("no"!=readConfStr("projPlot", "no", _runSec)) {
05458 cout<<endl<<" Do projPlot Action"<<endl<<endl;
05459
05460 initParams("ProjPlot", fullParams, fullParamsWOI, "yes", "mlFit", "no",
05461 kTRUE, postPdfFloatSet, preMLFixSet, preMLFloatSet);
05462
05463 cout << "Floating all parameters so that projection plots display correctly." << endl;
05464 fullParams.setAttribAll("Constant",kFALSE);
05465
05466 RooDataSet *projData=
05467 _datasets->getData(readConfStrCnA("projPlotData", "notSet"));
05468 if(!projData) {
05469 cout<<endl<<" Please specify in section ["<<_runSec<<"]"<<endl
05470 <<" the dataset you want to project with config "<<endl<<endl
05471 <<" projPlotData = <datasetName> // and (see online doc)"<<endl
05472 <<" projPlotData_<obsName> = <datasetName>"<<endl
05473 <<endl<<" and ru-run your job"<<endl;
05474 exit(-1);
05475 }
05476 chkBlind(projData->GetName());
05477
05478 TList *plotList=new TList();
05479
05480 getSnB();
05481
05482 TString varNames=readConfStr("projVars", "", _runSec);
05483 if (""==varNames) varNames=readConfStr("projVar","", _runSec);
05484 rarStrParser varNamesParser=varNames;
05485 for (Int_t i=0; i<varNamesParser.nArgs(); i++) {
05486
05487 RooRealVar *theVar=(RooRealVar*)_fullObs->find(varNamesParser[i]);
05488 if(!theVar) {
05489 cout<<"Can not find obs named "<<varNamesParser[i]<<endl;
05490 continue;
05491 }
05492 doProjPlot(projData, theVar, *plotList);
05493 }
05494
05495 doLLRPlot(projData, *plotList);
05496
05497 TString projPlotFile=readConfStr("projPlotFile", "default", _runSec);
05498 projPlotFile=getRootFileName("projPlot", projPlotFile);
05499 TFile plotFile(projPlotFile,"recreate");
05500 plotList->Print();
05501 plotList->Write();
05502 plotFile.Close();
05503 cout<<"Writing out proj plots to "<<projPlotFile<<endl;
05504 }
05505
05506
05507 if ("no"!=readConfStr("scanPlot", "no", _runSec)) {
05508 cout<<endl<<" Do scanPlot Action"<<endl<<endl;
05509
05510 initParams("ScanPlot", fullParams, fullParamsWOI, "yes", "mlFit", "no",
05511 kTRUE, postPdfFloatSet, preMLFixSet, preMLFloatSet);
05512
05513 TList *plotList=new TList();
05514 doScanPlot(*plotList);
05515 TString scanPlotFile=readConfStr("scanPlotFile", "default", _runSec);
05516 scanPlotFile=getRootFileName("scanPlot", scanPlotFile);
05517 TFile plotFile(scanPlotFile,"recreate");
05518 plotList->Print();
05519 plotList->Write();
05520 plotFile.Close();
05521 cout<<endl<<"Writing out scan plots (dataset) to "<<scanPlotFile<<endl;
05522 }
05523
05524
05525 if ("no"!=readConfStr("contourPlot", "no", _runSec)) {
05526 cout<<endl<<" Do contourPlot Action"<<endl<<endl;
05527
05528 initParams("ContourPlot", fullParams, fullParamsWOI, "yes", "mlFit", "no",
05529 kTRUE, postPdfFloatSet, preMLFixSet, preMLFloatSet);
05530
05531 TList *plotList=new TList();
05532 doContourPlot(*plotList);
05533 TString contourPlotFile=readConfStr("contourPlotFile", "default", _runSec);
05534 contourPlotFile=getRootFileName("contPlot", contourPlotFile);
05535 TFile plotFile(contourPlotFile,"recreate");
05536 plotList->Print();
05537 plotList->Write();
05538 plotFile.Close();
05539 cout<<"Writing out contour plots to "<<contourPlotFile<<endl;
05540 }
05541
05542
05543 if ("no"!=readConfStr("sPlot", "no", _runSec)) {
05544 cout<<endl<<" Do sPlot Action"<<endl<<endl;
05545
05546 initParams("SPlot", fullParams, fullParamsWOI, "yes", "mlFit", "no",
05547 kTRUE, postPdfFloatSet, preMLFixSet, preMLFloatSet);
05548
05549 TList *plotList=new TList();
05550
05551 TString varNames=readConfStr("sPlotVars", "", _runSec);
05552 if (""==varNames) varNames=readConfStr("sPlotVar","", _runSec);
05553 rarStrParser varNamesParser=varNames;
05554 for (Int_t i=0; i<varNamesParser.nArgs(); i++) {
05555
05556 RooRealVar *theVar=(RooRealVar*)_fullObs->find(varNamesParser[i]);
05557 if(!theVar) {
05558 cout<<"Can not find obs named "<<varNamesParser[i]<<endl;
05559 continue;
05560 }
05561 doSPlot(theVar, *plotList);
05562 }
05563 TString sPlotFile=readConfStr("sPlotFile", "default", _runSec);
05564 sPlotFile=getRootFileName("sPlot", sPlotFile);
05565 TFile plotFile(sPlotFile, "recreate");
05566 plotList->Print();
05567 plotList->Write();
05568 plotFile.Close();
05569 cout<<"Writing out sPlots to "<<sPlotFile<<endl;
05570 }
05571
05572
05573 if ("no"!=readConfStr("combinePlot", "no", _runSec)) {
05574 cout<<endl<<" Do combinePlot Action"<<endl<<endl;
05575
05576
05577
05578
05579 TList *plotList=new TList();
05580 doCombinePlot(*plotList);
05581 TString combinePlotFile=readConfStr("combinePlotFile", "default", _runSec);
05582 combinePlotFile=getRootFileName("combinePlot", combinePlotFile);
05583 TFile plotFile(combinePlotFile,"recreate");
05584 plotList->Print();
05585 plotList->Write();
05586 plotFile.Close();
05587 cout<<"Writing out combine plots to "<<combinePlotFile<<endl;
05588 }
05589 }
