rarBasePdf.cc

Go to the documentation of this file.
00001 /*****************************************************************************
00002 * Project: BaBar detector at the SLAC PEP-II B-factory
00003 * Package: RooRarFit
00004  *    File: $Id: rarBasePdf.cc,v 1.64 2013/09/11 11:58:21 fwilson Exp $
00005  * Authors: Lei Zhang
00006  * History:
00007  * 
00008  * Copyright (C) 2005-2012, University of California, Riverside
00009  *****************************************************************************/
00010 
00011 // -- CLASS DESCRIPTION [RooRarFit] --
00012 // This class provides a base Pdf class for RooRarFit
00014 //
00015 // BEGIN_HTML
00016 // This class provides a base Pdf class for RooRarFit
00017 // END_HTML
00018 //
00019 
00020 #include "RooRarFit/rarVersion.hh"
00021 
00022 #include "Riostream.h"
00023 #include <fstream>
00024 #include <sstream>
00025 using namespace std;
00026 
00027 #include "TArrayI.h"
00028 #include "TPaveText.h"
00029 #include "TLatex.h"
00030 #include "TAxis.h"
00031 #include "TH1F.h"
00032 
00033 #include "RooFitCore/RooAbsPdf.hh"
00034 #include "RooFitCore/RooAddPdf.hh"
00035 #include "RooFitCore/RooArgList.hh"
00036 #include "RooFitCore/RooArgSet.hh"
00037 #include "RooFitCore/RooDataSet.hh"
00038 #include "RooFitCore/RooFitResult.hh"
00039 #include "RooFitCore/RooHist.hh"
00040 #include "RooFitCore/RooPlot.hh"
00041 #include "RooFitCore/RooProdPdf.hh"
00042 #include "RooFitCore/RooRandom.hh"
00043 #include "RooFitCore/RooRealVar.hh"
00044 #include "RooFitCore/RooStringVar.hh"
00045 #include "RooFitCore/RooSuperCategory.hh"
00046 #include "RooFitCore/RooGlobalFunc.hh"
00047 using namespace RooFit;
00048 
00049 #include "RooRarFit/rarBasePdf.hh"
00050 #include "RooRarFit/rarMLFitter.hh"
00051 
00052 ClassImp(rarBasePdf)
00053   ;
00054 
00055 RooGenericPdf rarBasePdf::_dummyPdf("dummyPdf","dummy Pdf","1",RooArgList());
00056 RooConstVar rarBasePdf::_dummyExpEvt("dummyExpEvt", "dummyExpEvt", 0);
00057 RooExtendPdf rarBasePdf::_dummyExtPdf("dummyExtPdf","dummy extend Pdf",
00058                              rarBasePdf::_dummyPdf, rarBasePdf::_dummyExpEvt);
00059 RooCategory rarBasePdf:: _compCat("compCat", "compCat");
00060 rarMLFitter *rarBasePdf::_theFitter(0);
00061 // plot colors
00062 Int_t rarBasePdf::_rarColors[NCOLORS]={
00063   8, // light green
00064   6, // magenta
00065   9, // light blue
00066   49, // dark brown
00067   42, // brown
00068   //11, // light brown
00069   30, // dark green
00070   13, // dark gray
00071   7, // cyan
00072   2, // red
00073   3, // green
00074   5, // yellow
00075 };
00076 TString rarBasePdf::_fracNames;
00077 
00081 rarBasePdf::rarBasePdf()
00082   : rarConfig(),
00083     _pdfType("notSet"), _datasets(0), _theData(0),
00084     _thePdf(0), _theProtGen(0),
00085     _thisSimPdf(0), _thisSimPdfWOP(0), _theSimPdf(0),
00086     _nxPdf(0), _myDummyPdf(0), _controlStr("")
00087 {
00088   init();
00089 }
00090 
00103 rarBasePdf::rarBasePdf(const char *configFile, const char *configSec,
00104                        const char *configStr,
00105                        rarDatasets *theDatasets, RooDataSet *theData,
00106                        const char *name, const char *title)
00107   : rarConfig(configFile, configSec, configStr, name, title),
00108     _pdfType("notSet"), _datasets(theDatasets), _theData(theData),
00109     _thePdf(0), _theProtGen(0),
00110     _thisSimPdf(0), _thisSimPdfWOP(0), _theSimPdf(0),
00111     _nxPdf(0), _myDummyPdf(0), _controlStr("")
00112 {
00113   init();
00114 }
00115 
00116 rarBasePdf::~rarBasePdf()
00117 {
00118   // _dataSets.Delete();
00119 }
00120 
00144 void rarBasePdf::init()
00145 {
00146   // set pdf type
00147   rarStrParser configStrParser=_configStr;
00148   configStrParser.Remove(); // name has been set
00149   _pdfType=configStrParser[0];
00150   configStrParser.Remove();
00151   // if this is the main fitter? if none set yet, it should be
00152   if (("MLFitter"==_pdfType)&&(!getFitter())) setFitter((rarMLFitter*)this);
00153   
00154   // set paramSec
00155   TString paramSec=readConfStr("paramSec_"+_pdfType, "notSet");
00156   if ("notSet"==paramSec) {
00157     if ("notSet"==(paramSec=readConfStr("paramSec", "notSet"))) {
00158       paramSec=_configSec;
00159     }
00160   }
00161   setVarSec(paramSec);
00162   cout<<"init of rarBasePdf:"<<endl
00163       <<" pdf Type: "<<_pdfType<<endl
00164       <<" pdfTitle: "<<GetTitle()<<endl
00165       <<" paramSec: "<<getVarSec()<<endl;
00166   
00167   // set control bits
00168   setControlBit("PdfFit", "pdfFit");
00169   setControlBit("PdfPlot", "pdfPlot");
00170   setControlBit("ParamsOnPlot", "paramsOnPlot");
00171   setControlBit("Chi2OnPlot", "chi2OnPlot");
00172   setControlBit("FirstFitOnly", "firstFitOnly");
00173   // default no
00174   setControlBit("noCompsOnPlot", "compsOnPlot");
00175   setControlBit("noCompsDataOnPlot", "compsDataOnPlot");
00176   // for SimPdf
00177   setControlBit("noSimFit", "simultaneousFit", getMasterSec());
00178   
00179   //setControlBit("noUseGenOnly", "useGeneratorOnly", getMasterSec());
00180   
00181   // get full obs
00182   _fullObs=_datasets->getFullFObs();
00183   // get the default dataset
00184   setFitData(_theData);
00185   
00186   // creates extra params if needed
00187   // The extra params created may or my not be directly used by this Pdf
00188   createAbsVars("xtraParams", &_xParams);
00189   
00190   // creates extra pdfs if needed
00191   // It creates extra PDFs associated with this RooRarFitPdf.
00192   // The extra pdfs themselves might not be direct components
00193   // of the (final) PDFs, but their parameters, components, etc.,
00194   // can be part of the (final) PDFs.
00195   createPdfs("xtraPdfs", &_xPdfList);
00196   // creates extra generators if needed
00197   createPdfs("xtraGenerators", &_xPdfList, &_protGenPdfs);
00198   
00199   // add protDataVars
00200   addProtVars();
00201 
00202 }
00203 
00209 void rarBasePdf::addProtVars()
00210 {
00211   addProtVars("conditionalObs", _conditionalObs);
00212   addProtVars("protDataVars", _protDataVars);
00213   _protDataVars.add(_conditionalObs, kTRUE);
00214   // add _conditionalObs for pdfFit
00215   _condObss.add(_conditionalObs);
00216 }
00217 
00223 void rarBasePdf::addProtVars(TString configName, RooArgSet &protVars)
00224 {
00225   // get protDataVars defined in config file
00226   RooArgSet allProtVars(getArgSet(configName, kTRUE));
00227   allProtVars.add(getArgSet(readConfStr(configName,"", _runSec),kFALSE));
00228   if (allProtVars.getSize()>0) {
00229     //cout<<"protDataVars (conditionalObs) defined:"<<endl;
00230     //allProtVars.Print("v");
00231     RooArgList allProtVarList(allProtVars);
00232     for (Int_t i=0; i<allProtVarList.getSize(); i++) {
00233       TString varName=allProtVarList[i].GetName();
00234       RooAbsArg *theVar=_fullObs->find(varName);
00235       if (!theVar) {
00236         cout<<" Can not find obs named "<<varName<<endl;
00237       }
00238       protVars.add(*theVar);
00239     }
00240   }
00241 }
00242 
00248 void rarBasePdf::setCondObss(RooArgSet condObsSet)
00249 {
00250   if (condObsSet.getSize()<=0) return;
00251   if (!_thePdf) return;
00252   if (!_theData) return;
00253   // first get all observables
00254   _condObss.add(*(_thePdf->getObservables(_theData)));
00255   // remove CondObss
00256   _condObss.remove(condObsSet);
00257   cout<<" conditionalObs for pdffit"<<endl;
00258   _condObss.Print("v");
00259 }
00260 
00269 void rarBasePdf::setFitData(RooDataSet *theData)
00270 {
00271   _theData=theData;
00272   // check if fitData is defined in param section
00273   TString dsName=readConfStr("fitData", "notSet", getVarSec());
00274   if ("notSet"!=dsName) {
00275     RooDataSet *myData=_datasets->getData(dsName);
00276     if (myData) _theData=myData;
00277     else {
00278       cout<<" W A R N I N G !"<<endl
00279           <<" Dataset named "<<dsName<<" can not be found"<<endl
00280           <<" fitData = "<<dsName<<endl
00281           <<" defined in section "<<getVarSec()<<endl;
00282       _theData=0;
00283     }
00284   }
00285   
00286   if (_theData) {
00287     cout<<" default Dataset for "<<GetName()<<": "
00288         <<_theData->GetName()<<endl;
00289   } else {
00290     cout<<" no default dataset for "<<GetName()<<endl;
00291   }
00292   return;
00293 }
00294 
00300 RooArgSet rarBasePdf::getParams()
00301 {
00302   RooArgSet theSet(_params);
00303   for (Int_t i=0; i<_nxPdf; i++) {
00304     rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
00305     theSet.add(thePdf->getParams());
00306   }
00307   theSet.add(_xParams);
00308   
00309   return theSet;
00310 }
00311 
00317 void rarBasePdf::addToParams(RooRealVar *theVar)
00318 {
00319   _params.add(*theVar);
00320 }
00321 
00327 void rarBasePdf::addToObs(RooRealVar *theVar)
00328 {
00329   _obsSet.add(*theVar);
00330 }
00331 
00343 RooArgList *rarBasePdf::getFormulaArgs(rarStrParser fStrParser)
00344 {
00345   RooArgList *depList=new RooArgList;
00346   Int_t nArgs=fStrParser.nArgs();
00347   if (nArgs<=0) return depList;
00348   for (Int_t i=0; i<nArgs; i++) {
00349     RooAbsArg *theDep(0);
00350     if (theDep=_obsSet.find(fStrParser[i]));
00351     else if (theDep=_params.find(fStrParser[i]));
00352     else if (theDep=_fullObs->find(fStrParser[i])) {
00353       _obsSet.add(*theDep);
00354     } else { // create it
00355       if (isNumber(fStrParser[i])) break;
00356       theDep=createAbsReal(fStrParser[i], fStrParser[i]);
00357     }
00358     depList->add(*theDep);
00359   }
00360   
00361   return depList;
00362 }
00363 
00373 Double_t rarBasePdf::getFormulaVal(TString varStr)
00374 {
00375   Double_t retVal(0);
00376   rarStrParser varStrParser=varStr;
00377   TString formula=varStrParser[0]; // get the formula
00378   varStrParser.Remove();
00379   // we need to prserve _params RooArgSet
00380   RooArgSet params(_params);
00381   RooFormulaVar myVar("myVar", formula, *getFormulaArgs(varStrParser));
00382   retVal=myVar.getVal();
00383   // restore original _params
00384   _params.removeAll();
00385   _params.add(params);
00386   return retVal;
00387 }
00388 
00398 rarBasePdf *rarBasePdf::createPdfs(TString Comps, TList *pdfList,
00399                                    RooAbsCollection *PDFs, TString secName)
00400 {
00401   // secName
00402   if (""==secName) secName=getVarSec();
00403   // get number of components
00404   rarStrParser compStrParser=readConfStr(Comps, "", secName);
00405   Int_t nComp=compStrParser.nArgs();
00406   // read in comp info
00407   RooArgSet pdfCompSet("PDF Comp Set");
00408   for (Int_t i=0; i<nComp; i++) {
00409     TString configStr=
00410       readConfStr("configStr", "notSet", compStrParser[i]+" Config");
00411     RooStringVar *pdf=
00412       new RooStringVar(compStrParser[i], compStrParser[i], configStr);
00413     pdfCompSet.addOwned(*pdf);
00414   }
00415   //pdfCompSet.Print("v");
00416   // list configed pdfs
00417   if (nComp>0) {
00418     cout <<endl<<"Pdfs defined with config \""<<Comps<<"\""
00419          <<" in config file for "<<GetName()<<":"<<endl;
00420   }
00421   for (Int_t i=0; i<nComp; i++) {
00422     TString pdfStr=(RooStringVar&)pdfCompSet[compStrParser[i]];
00423     cout<<Form(" #%02d ",i)<<compStrParser[i]<<" "<<pdfStr<<endl;
00424   }
00425   cout<<endl;
00426   // now build individual pdf
00427   rarBasePdf *thePdf(0);
00428   for (Int_t i=0; i<nComp; i++) {
00429     TString pdfStr=(RooStringVar&)pdfCompSet[compStrParser[i]];
00430     thePdf=createPdf(compStrParser[i]+" "+pdfStr);
00431     if (PDFs) PDFs->add(*thePdf->getPdf());
00432     if (pdfList) pdfList->Add(thePdf);
00433   }
00434   
00435   return thePdf;
00436 }
00437 
00450 void rarBasePdf::preAction()
00451 {
00452   // first for xtraPdf
00453   _nxPdf=_xPdfList.GetSize();
00454   for (Int_t i=0; i<_nxPdf; i++) {
00455     rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
00456     thePdf->preAction();
00457   }
00458   cout<<" In rarBasePdf preAction for "<<GetName()<<endl;
00459   // first get full obs including those from sub pdfs
00460   _fObsSet.add(*_thePdf->getDependents(*_fullObs));
00461   // create _myDummyPdf
00462   if (!_myDummyPdf) {
00463     RooArgSet myFullVars(*_thePdf->getParameters(_theData));
00464     _myDummyPdf=new RooGenericPdf
00465       (Form("myDummyPdf_%s", GetName()), "1", myFullVars);
00466   }
00467   // build simPdf for this Pdf
00468   _thisSimPdf=(RooSimultaneous*)getSimPdf(_theSimPdf);
00469   _thisSimPdfWOP=_thisSimPdf;
00470 }
00471 
00478 RooAbsPdf *rarBasePdf::getSimPdf(RooSimultaneous *simPdf, RooAbsPdf *srcPdf)
00479 {
00480   if (!simPdf) return _thisSimPdf;
00481   RooSimultaneous *theSimPdf(0);
00482   if (TString(simPdf->ClassName())!="RooSimultaneous") return theSimPdf;
00483   if (!srcPdf) srcPdf=_thePdf;
00484   // simPdf for master pdf is itself
00485   if ((simPdf==_theSimPdf)&&(srcPdf==_theSimPdf)) return _theSimPdf;
00486   // build simPdf for the pdf based on full SimPdf
00487   RooAbsCategoryLValue *simCats=(RooAbsCategoryLValue *)(&simPdf->indexCat());
00488   //simCats->Print("v");
00489   // get pdf name
00490   theSimPdf=new RooSimultaneous(Form("simPdf_%s", srcPdf->GetName()),
00491                                 Form("simPdf for %s", srcPdf->GetTitle()),
00492                                 *simCats);
00493   Bool_t isFound=kFALSE;
00494   Int_t nCats=simCats->numTypes();
00495   for(Int_t i=0; i<nCats; i++) {
00496     RooCatType *catType=(RooCatType*)simCats->lookupType(i);
00497     TString catName=catType->GetName();
00498     TString catName2=catName;
00499     // do we have physCat?
00500     TString physCat=getFitter()->getPhysCat();
00501     if ((""!=physCat)&&(catName2.Contains(";"))) {
00502       TString physCatType=catName2(1,catName2.First(';')-1);
00503       catName2.ReplaceAll(physCatType+";", "");
00504       catName2.ReplaceAll("}", ";"+physCatType+"}");
00505     }
00506     // check if we can find cloned pdf
00507     RooAbsPdf *thisCatFullPdf=simPdf->getPdf(catName.Data());
00508     if (!thisCatFullPdf) {
00509       theSimPdf->addPdf(_dummyPdf, catName);
00510     } else {
00511       //thisCatFullPdf->Print("v");
00512       RooArgSet *thisCatCompSet=thisCatFullPdf->getComponents();
00513       RooAbsPdf *thisCatPdf=(RooAbsPdf*)thisCatCompSet->
00514         find(Form("%s_%s", srcPdf->GetName(), catName.Data()));
00515       if (!thisCatPdf) thisCatPdf=(RooAbsPdf*)thisCatCompSet->
00516         find(Form("%s_%s", srcPdf->GetName(), catName2.Data()));
00517       if (thisCatPdf) {
00518         theSimPdf->addPdf(*thisCatPdf, catName);
00519         isFound=kTRUE;
00520       } else {
00521         theSimPdf->addPdf(_dummyPdf, catName);
00522       }
00523     }
00524   }
00525   // do we really find sim pdf for this pdf
00526   if (!isFound) {
00527     // this Pdf does not have SimPdf counter-part
00528     //cout<<" No need to build simPdf for "<<GetName()<<endl;
00529     delete theSimPdf;
00530     theSimPdf=0;
00531   }
00532   if (theSimPdf) {
00533     cout<<" simPdf created for "<<srcPdf->GetName()<<endl;
00534     theSimPdf->Print();
00535     //theSimPdf->Print("v");
00536   }
00537   return theSimPdf;
00538 }
00539 
00555 RooArgSet rarBasePdf::getArgSet(TString paramNames, Bool_t useRead,
00556                                 RooArgSet *fullSet)
00557 {
00558   RooArgSet retVal;
00559   rarStrParser paramNameStrParser="";
00560   if (!useRead) {
00561     paramNameStrParser=paramNames;
00562   } else {
00563     paramNameStrParser=readConfStr(paramNames, "", getVarSec());
00564   }
00565   Int_t nParams=paramNameStrParser.nArgs();
00566   RooAbsArg *theParam(0);
00567   rarBasePdf *thePdf(0);
00568   for (Int_t i=0; i<nParams; i++) {
00569     TString paramName=paramNameStrParser[i];
00570     // if fullSet have the name
00571     if (fullSet) {
00572       RooArgSet *inFS=(RooArgSet*)fullSet->selectByName(paramName);
00573       retVal.add(*inFS);
00574       delete inFS;
00575     }
00576     // check if the paramName is the pdf Name
00577     if(paramName==GetName()) {
00578       retVal.add(_params);
00579       // have any split param?
00580       if (fullSet) {
00581         RooArgList paramList(_params);
00582         for (Int_t j=0; j<paramList.getSize(); j++) {
00583           RooArgSet *inFS=(RooArgSet*)
00584             fullSet->selectByName(Form("%s_*",paramList[j].GetName()));
00585           retVal.add(*inFS);
00586           delete inFS;
00587         }
00588       }
00589     } else if (theParam=getAbsVar(paramName)) {
00590       //RooAbsArg *theParam=_params.find(paramName);
00591       retVal.add(*theParam);
00592       if (fullSet) {
00593         RooArgSet *inFS=(RooArgSet*)
00594           fullSet->selectByName(Form("%s_*",theParam->GetName()));
00595         retVal.add(*inFS);
00596         delete inFS;
00597       }
00598     } else if (thePdf=(rarBasePdf*)(_rarPdfs.FindObject(paramName))) {
00599       retVal.add(thePdf->getArgSet(thePdf->GetName(), kFALSE, fullSet));
00600     }
00601   }
00602   // now for xtraPdf
00603   for (Int_t i=0; i<_nxPdf; i++) {
00604     rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
00605     retVal.add(thePdf->getArgSet(paramNames, useRead, fullSet));
00606   }
00607   return retVal;
00608 }
00609 
00614 Bool_t rarBasePdf::matchCatType(RooCatType *catTypeN, RooCatType *catTypeO)
00615 {
00616   Bool_t retVal(kTRUE);
00617   TString typeNameO=catTypeO->GetName();
00618   typeNameO.ReplaceAll("{", "");
00619   typeNameO.ReplaceAll("}", "");
00620   typeNameO.ReplaceAll(";", " ");
00621   rarStrParser typeParserO=typeNameO;
00622   
00623   TString typeNameN=catTypeN->GetName();
00624   typeNameN.ReplaceAll("{", "");
00625   typeNameN.ReplaceAll("}", "");
00626   typeNameN.ReplaceAll(";", " ");
00627   rarStrParser typeParserN=typeNameN;
00628   for (Int_t i=0; i<typeParserN.nArgs(); i++) {
00629     if (!typeParserO.Have(typeParserN[i])) {
00630       retVal=kFALSE;
00631       break;
00632     }
00633   }
00634   
00635   return retVal;
00636 }
00637 
00642 void rarBasePdf::saveFracName(TString fracName)
00643 {
00644   if (isFracName(fracName)) return;
00645   _fracNames+=" "+fracName;
00646   //cout<<_fracNames<<endl;
00647 }
00648 
00654 Bool_t rarBasePdf::isFracName(TString fracName)
00655 {
00656   Bool_t retVal(kFALSE);
00657   rarStrParser fracNamesParser=_fracNames;
00658   if (fracNamesParser.Have(fracName)) retVal=kTRUE;
00659   
00660   return retVal;
00661 }
00662 
00667 void rarBasePdf::attachDataSet(const RooAbsData &data)
00668 {
00669   RooAbsPdf *thePdf(0);
00670   if (_thePdf) thePdf=_thePdf;
00671   if (_thisSimPdf) thePdf=_thisSimPdf;
00672   if (!thePdf) return;
00673   thePdf->attachDataSet(data);
00674 }
00675 
00678 Bool_t rarBasePdf::isNegativeValue()
00679 {
00680   Bool_t isNegative(kFALSE);
00681   RooAbsPdf *thePdf(0);
00682   if (_thePdf) thePdf=_thePdf;
00683   if (_thisSimPdf) thePdf=_thisSimPdf;
00684   if (!thePdf) return isNegative;
00685   if (thePdf->getVal()<0) {
00686     isNegative=kTRUE;
00687     cout<<" Pdf "<<thePdf->GetName()<<" has negative value"<<endl;
00688     thePdf->Print();
00689     thePdf->Print("v");
00690     return isNegative;
00691   }
00692   return isNegative;
00693 }
00694 
00695 
00707 RooBinning *rarBasePdf::getRange(RooRealVar *theVar, TString rPrefix,
00708                                  Double_t &min, Double_t &max,
00709                                  const Char_t *sec, Int_t *nBins)
00710 {
00711   RooBinning *abins(0);
00712   if (!sec) sec=getVarSec();
00713   min=theVar->getMin();
00714   max=theVar->getMax();
00715   rarStrParser rangeParser=
00716     readConfStr(rPrefix+theVar->GetName(),"", sec);
00717   if (rangeParser.nArgs()<=0) { // try runSec then
00718     rangeParser=readConfStr(rPrefix+theVar->GetName(),"", _runSec);
00719   }
00720   if (rangeParser.nArgs()>0) {
00721     min=atof(rangeParser[0]);
00722     rangeParser.Remove();
00723     if (min<theVar->getMin()) min=theVar->getMin();
00724   }
00725   if (rangeParser.nArgs()>0) {
00726     max=atof(rangeParser[rangeParser.nArgs()-1]);
00727     rangeParser.Remove(rangeParser.nArgs()-1);
00728     if (max>theVar->getMax()) max=theVar->getMax();
00729   }
00730   if (max<min) {
00731     Double_t v=max;
00732     max=min;
00733     min=v;
00734   }
00735   if (nBins) {
00736     abins=new RooBinning(*nBins, min, max);
00737     if (rangeParser.nArgs()>0) {
00738       delete abins;
00739       abins=new RooBinning(min, max);
00740       while (rangeParser.nArgs()>0) {
00741         Double_t boundary=atof(rangeParser[0]);
00742         rangeParser.Remove();
00743         if ((boundary<min) || (boundary>max)) continue;
00744         abins->addBoundary(boundary);
00745       }
00746     }
00747   }
00748   return abins;
00749 }
00750 
00757 RooArgSet rarBasePdf::getCorrCoeffs()
00758 {
00759   RooArgSet theSet(_corrCoeffs);
00760   for (Int_t i=0; i<_nxPdf; i++) {
00761     rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
00762     theSet.add(thePdf->getCorrCoeffs());
00763   }
00764   
00765   return theSet;
00766 }
00767 
00775 Double_t rarBasePdf::getCorrCoeff(const TString pn1, const TString pn2)
00776 {
00777   // trivial case, but need to be considered specially
00778   if (pn1==pn2) return 1;
00779   // do we have it in argset?
00780   RooArgSet theCorrSet(getCorrCoeffs());
00781   return theCorrSet.getRealValue(getCorrCoefName(pn1, pn2), 0);
00782 }
00783 
00788 void rarBasePdf::saveCorrCoeffs(RooFitResult *fr)
00789 {
00790   // add correlation matrix
00791   Int_t nFloat=fr->floatParsFinal().getSize();
00792   for (Int_t iF=0; iF<nFloat; iF++) {
00793     RooAbsArg &iArg=fr->floatParsFinal().operator[](iF);
00794     for(Int_t jF=0; jF<iF; jF++) {
00795       RooAbsArg &jArg=fr->floatParsFinal().operator[](jF);
00796       TString corrName=getCorrCoefName(iArg.GetName(), jArg.GetName());
00797       saveCorrCoeff(corrName, fr->correlation(iArg, jArg));
00798     }
00799   }
00800 }
00801 
00809 Bool_t rarBasePdf::saveCorrCoeff(TString corrName, Double_t corrCoef,
00810                                  Bool_t saveTrivial)
00811 {
00812   Bool_t retVal(kFALSE);
00813   if (fabs(corrCoef)<1e-4) {
00814     if (!saveTrivial)
00815       return retVal;
00816   }
00817   if (fabs(corrCoef)>1) return retVal;
00818   if (_corrCoeffs.setRealValue(corrName, corrCoef)) {
00819     // create corr coeff RRV
00820     RooRealVar *theCorrCoefVar=new RooRealVar(corrName, corrName, corrCoef);
00821     _corrCoeffs.add(*theCorrCoefVar);
00822   }
00823   retVal=kTRUE;
00824   
00825   return retVal;
00826 }
00827 
00832 TString rarBasePdf::getCorrCoefName(const TString pn1, const TString pn2) const
00833 {
00834   // construct corrcoef var name
00835   TString corrName=Form("corrCoef_%s_%s", pn1.Data(), pn2.Data());
00836   if (pn1.CompareTo(pn2)>0)
00837     corrName=Form("corrCoef_%s_%s", pn2.Data(), pn1.Data());
00838   return corrName;
00839 }
00840 
00845 void rarBasePdf::doXPdfFit(TString pdfList)
00846 {
00847   for (Int_t i=0; i<_nxPdf; i++) {
00848     rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
00849     thePdf->doPdfFit(pdfList);
00850   }
00851 }
00852 
00866 void rarBasePdf::doPdfFit(TString pdfList)
00867 {
00868   // first for xtraPfs
00869   doXPdfFit(pdfList);
00870   
00871   cout<<endl<<" In rarBasePdf doPdfFit for "<<GetName()<<endl;
00872   
00873   if (!getControlBit("PdfFit")) return;
00874   TString pdfFitStr=readConfStr("pdfFit","yes",getVarSec());
00875   // is fit been done?
00876   if (getControlBit("FirstFitOnly")) {
00877     if (getControlBit("FirstFitDone")) return;
00878   }
00879   // are we in the list for plotting
00880   rarStrParser pdfListParser=pdfList;
00881   if (pdfListParser.nArgs()>0)
00882     if (!pdfListParser.Have(GetName())) return;
00883   
00884   // first check if it has its pdf and data
00885   assert(_thePdf);
00886   assert(_theData);
00887   // display dataset to fit
00888   cout<<" pdfFit to dataset: "<<_theData->GetName()<<endl;
00889   // fitoption
00890   TString fitOption="hrq";
00891   //if (_thePdf->isExtended()) fitOption="emhr";
00892   // convert to newer fitTo format for steering options
00893   Bool_t pdfExtended = fitOption.Contains("e");
00894   Bool_t pdfMinos    = fitOption.Contains("m");
00895   Bool_t pdfHesse    = fitOption.Contains("h");
00896   Bool_t pdfVerbose  = !fitOption.Contains("q");
00897   Bool_t pdfSave     = fitOption.Contains("r"); // return results
00898  
00899   // get number of cpus option
00900   Int_t pdfFitNumCPU=atoi(readConfStr("useNumCPU", "1", getMasterSec()));
00901 
00902   // save obs
00903   string obsSaveStr;
00904   writeToStr(_fObsSet, obsSaveStr);
00905   // set fit ranges
00906   RooArgList fObsList(_fObsSet);
00907   for (Int_t i=0; i<fObsList.getSize(); i++) {
00908     // get obs
00909     RooRealVar *theObs=dynamic_cast<RooRealVar*>(fObsList.at(i));
00910     if (!theObs) continue;
00911     Double_t fitMin=theObs->getMin();
00912     Double_t fitMax=theObs->getMax();
00913     getRange(theObs, "fitRange_", fitMin, fitMax);
00914     theObs->setRange(fitMin, fitMax);
00915   }
00916   
00917   // get all params for this pdf and simPdf
00918   RooArgSet thisFParams(*_thePdf->getParameters(_theData));
00919   if (_thisSimPdf) thisFParams.add(*_thisSimPdf->getParameters(_theData));
00920   // first fix params specified with prePdfFix
00921   RooArgSet prePdfFixParamSet=
00922     getArgSet(readConfStr("prePdfFix","", getVarSec()), kFALSE, &thisFParams);
00923   // set those params to values specified
00924   rarStrParser prePdfFixParser=readConfStr("prePdfFix", "", getVarSec());
00925   while (prePdfFixParser.nArgs()>0) {
00926     TString varName=prePdfFixParser[0];
00927     prePdfFixParser.Remove();
00928     RooRealVar *theVar=(RooRealVar*)prePdfFixParamSet.find(varName);
00929     if (!theVar) continue;
00930     // check if the next is number
00931     TString varVal=prePdfFixParser[0];
00932     if (isNumber(varVal)) {
00933       prePdfFixParser.Remove();
00934       if (theVar->ClassName()==TString("RooRealVar"))
00935         theVar->setVal(atof(varVal));
00936     }
00937   }
00938   if (prePdfFixParamSet.getSize()>0) {
00939     cout<<" prePdfFix:"<<endl;
00940     prePdfFixParamSet.Print("v");
00941   }
00942   // then float params specified w/ prePdfFloat
00943   RooArgSet prePdfFloatParamSet=
00944     getArgSet(readConfStr("prePdfFloat","", getVarSec()),kFALSE,&thisFParams);
00945   // merge those params to prePdfFix for saving
00946   prePdfFixParamSet.add(prePdfFloatParamSet);
00947   string prePdfFixParamSetSSaver;
00948   writeToStr(prePdfFixParamSet, prePdfFixParamSetSSaver);
00949   // fix or float
00950   prePdfFixParamSet.setAttribAll("Constant"); 
00951   prePdfFloatParamSet.setAttribAll("Constant", kFALSE);
00952   
00953   // fit to default dataset
00954   if ("simFitOnly"!=pdfFitStr) {
00955     cout<<endl<<" In rarBasePdf doPdfFit for "<<GetName()<< " Options: " << fitOption 
00956       << " using " << pdfFitNumCPU << " CPUs (if available)." << endl;
00957     RooFitResult *fitResult=
00958       _thePdf->fitTo(*_theData,ConditionalObservables(_condObss),
00959                      Save(pdfSave), Extended(pdfExtended), 
00960                      Verbose(pdfVerbose), Hesse(pdfHesse),
00961                      Minos(pdfMinos), NumCPU(pdfFitNumCPU));
00962 
00963     saveCorrCoeffs(fitResult);
00964   }
00965   // for simPdf
00966   if (_thisSimPdf&&pdfFitStr.Contains("simFit")) {
00967     cout<<"Should also have simPdf fit"<<endl;
00969     //string coeffSSaver;
00970     //writeToStr(_params, coeffSSaver);
00971     //_params.setAttribAll("Constant");
00973     // let's find if there is any protCat?
00974     RooArgSet protDataEVars(getFitter()->getProtDataEVars());
00975     RooAbsCategoryLValue *simCats=(RooAbsCategoryLValue *)
00976       (&_thisSimPdf->indexCat());
00977     RooSuperCategory *sCatO=dynamic_cast<RooSuperCategory*>(simCats);
00978     if (sCatO) {
00979       RooArgSet inputCats(sCatO->inputCatList());
00980       Int_t nCatsO=inputCats.getSize();
00981       inputCats.remove(protDataEVars, kFALSE, kTRUE);
00982       if ((inputCats.getSize()!=nCatsO)&&(inputCats.getSize()>0)) {
00983         RooSuperCategory *sCatN=new
00984           RooSuperCategory("sCatN", "sCatN", inputCats);
00985         _thisSimPdfWOP=new RooSimultaneous
00986           (Form("%sWOP", _thisSimPdfWOP->GetName()), "simPDFWOP", *sCatN);
00987         Int_t nTypesN=sCatN->numTypes();
00988         Int_t nTypesO=sCatO->numTypes();
00989         for(Int_t i=0; i<nTypesN; i++) {
00990           RooCatType *catTypeN=(RooCatType*)sCatN->lookupType(i);
00991           for (Int_t j=0; j<nTypesO; j++) {
00992             RooCatType *catTypeO=(RooCatType*)sCatO->lookupType(j);
00993             if (matchCatType(catTypeN, catTypeO)) {
00994               RooAbsPdf *typePdf=_thisSimPdf->getPdf(catTypeO->GetName());
00995               if (TString(typePdf->GetName())==_dummyPdf.GetName()) continue;
00996               _thisSimPdfWOP->addPdf(*typePdf, catTypeN->GetName());
00997               break;
00998             }
00999           }
01000         }
01001       }
01002     }
01003     //RooFitResult *fitResult=
01004     //  _thisSimPdfWOP->fitTo(*_theData,_condObss,fitOption);
01005     RooFitResult *fitResult= _thisSimPdfWOP->fitTo(*_theData,ConditionalObservables(_condObss),
01006                                                    Save(pdfSave), Extended(pdfExtended), 
01007                                                    Verbose(pdfVerbose), Hesse(pdfHesse),
01008                                                    Minos(pdfMinos), NumCPU(pdfFitNumCPU));
01009     saveCorrCoeffs(fitResult);
01011     //readFromStr(_params, coeffSSaver);
01012   }
01013   // after fit, save the params for later use, like plotting, etc
01014   RooArgSet *theParams=_thePdf->getParameters(_theData);
01015   writeToStr(*theParams, _afterFitSaverStr);
01016   // restore params in prePdfFixParamSet
01017   readFromStr(prePdfFixParamSet, prePdfFixParamSetSSaver);
01018   // restore obs saved
01019   readFromStr(_fObsSet, obsSaveStr);
01020   setControlBits("FirstFitDone");
01021 }
01022 
01031 RooPlot *rarBasePdf::doXPdfPlot(TList &plotList, TString pdfList)
01032 {
01033   RooPlot *frame(0);
01034   for (Int_t i=0; i<_nxPdf; i++) {
01035     rarBasePdf *thePdf=(rarBasePdf*)_xPdfList.At(i);
01036     frame=thePdf->doPdfPlot(plotList, pdfList);
01037   }
01038   return frame;
01039 }
01040 
01055 RooPlot *rarBasePdf::doPdfPlot(TList &plotList, TString pdfList)
01056 {
01057   RooPlot *frame(0);
01058   // first for xtraPdfs
01059   frame=doXPdfPlot(plotList, pdfList);
01060   
01061   cout<<endl<<" In rarBasePdf doPdfPlot for "<<GetName()<<endl;
01062   
01063   //if (!getControlBit("PdfFit")) return frame;
01064   TString pdfFitStr=readConfStr("pdfFit","yes",getVarSec());
01065   if (!getControlBit("PdfPlot")) return frame;
01066   if (getControlBit("PdfPlotDone")) return frame;
01067   // are we in the list for plotting
01068   rarStrParser pdfListParser=pdfList;
01069   if (pdfListParser.nArgs()>0)
01070     if (!pdfListParser.Have(GetName())) return frame;
01071   
01072   // check if it has its pdf and data
01073   if ((!_thePdf)||(!_theData)) return frame;
01074   
01075   // save obs
01076   string obsSaveStr;
01077   writeToStr(_fObsSet, obsSaveStr);
01078   // restore the fit value for plotting
01079   string coeffParamSSaver;
01080   RooArgSet *theParams=_thePdf->getParameters(_theData);
01081   writeToStr(*theParams, coeffParamSSaver);
01082   readFromStr(*theParams, _afterFitSaverStr);
01083   
01084   // plot for each obs
01085   RooArgList fObsList(_fObsSet);
01086   for (Int_t i=0; i<fObsList.getSize(); i++) {
01087     // get obs
01088     RooRealVar *theObs=dynamic_cast<RooRealVar*>(fObsList.at(i));
01089     if (!theObs) continue;
01090     TString nBinStr="0";
01091     // if obs is in _protDataVars or _condObss,
01092     // do not plot unless explicitly requested
01093     if (_protDataVars.find(theObs->GetName())) nBinStr="-1";
01094     if (_condObss.find(theObs->GetName())) nBinStr="-1";
01095     // get nBins
01096     Int_t nBins=atoi(readConfStr(Form("plotBins_%s", theObs->GetName()),
01097                                  nBinStr, getVarSec()));
01098     if (nBins<0) continue;
01099     if (0==nBins) nBins=theObs->getBins();
01100     // plot range
01101     Double_t plotMin=theObs->getMin();
01102     Double_t plotMax=theObs->getMax();
01103     getRange(theObs, "plotRange_", plotMin, plotMax);
01104     // normalization CmdArg
01105     // see post on HN from Wouter for another way
01106     // to handle fit and plot for sub-ranges
01107     // using named ranges, like:
01108     // mes.setRange("signal", 5.27, 5.29);
01109     // de.setRange("signal",-0.1,0.1) ;
01110     // data->plotOn(frame2,CutRange("signal")) ;
01111     // model.plotOn(frame2,ProjectionRange("signal")) ;
01112     RooCmdArg plotCutNormCmd;
01113     TString rangeCuts="1";
01114     if (plotMin>theObs->getMin())
01115       rangeCuts+=Form("&&(%s>%f)",theObs->GetName(), plotMin);
01116     if (plotMax<theObs->getMax())
01117       rangeCuts+=Form("&&(%s<%f)",theObs->GetName(), plotMax);
01118     if ("1"!=rangeCuts) {
01119       Double_t d=_theData->sumEntries();
01120       Double_t n=_theData->sumEntries(rangeCuts);
01121       Double_t plotRangeNorm=n/d;
01122       cout<<" Based on total evts "<<d<<" and fitted evts "<<n<<endl
01123           <<" norm. scale factor for plotting: "<<plotRangeNorm<<endl
01124           <<endl;
01125       if (plotRangeNorm!=1) {
01126         cout<<" Using relative norm. factor "<<plotRangeNorm<<endl;
01127         plotCutNormCmd=Normalization(plotRangeNorm, RooAbsReal::Relative);
01128       }
01129     }
01130     theObs->setRange(plotMin, plotMax);
01131     // do we have projWData for plotting
01132     RooDataSet *projWData(0);
01133     TString projWDataStr=
01134       readConfStr(Form("projWData_%s", theObs->GetName()), "no", getVarSec());
01135     if ("yes"==projWDataStr) projWData=_theData;
01136     else if ("no"!=projWDataStr)
01137       projWData=getDatasets()->getData(projWDataStr);
01138     if (projWData) {
01139       cout<<" Get plot of "<<theObs->GetName()<<" for Pdf "<<_thePdf->GetName()
01140           <<" using projWData "<<projWData->GetName()<<endl;
01141     }
01142     // any plots for simPdf?
01143     if (_thisSimPdfWOP&&pdfFitStr.Contains("simFit")) {
01144       // all sim comps
01145       RooArgSet *simComps=(RooArgSet*)_thisSimPdfWOP->getComponents()->
01146         selectByName(Form("%s_*",_thePdf->GetName()));
01147       RooAbsCategoryLValue *simCats=(RooAbsCategoryLValue *)
01148         (&((RooSimultaneous*)_thisSimPdfWOP)->indexCat());
01149       RooArgSet inputCats;
01150       RooSuperCategory *simSuperCats=dynamic_cast<RooSuperCategory*>(simCats);
01151       if (simSuperCats) inputCats.add(simSuperCats->inputCatList());
01152       else inputCats.add(*simCats);
01153       TIterator *cIter=inputCats.createIterator();
01154       RooCategory *theCat(0);
01155       while(theCat=(RooCategory*)cIter->Next()) {
01156         TIterator *tIter=theCat->typeIterator();
01157         RooCatType *theType(0);
01158         while (theType=(RooCatType*)tIter->Next()) {
01159           // do we have simPdf comp for this cat?
01160           TString compStr="";
01161           RooArgSet compParams;
01162           TIterator *compIter=simComps->createIterator();
01163           RooAbsPdf *thePdf(0);
01164           while(thePdf=(RooAbsPdf*)compIter->Next()) {
01165             TString thePdfName=thePdf->GetName();
01166             thePdfName.Remove(0, strlen(_thePdf->GetName()));
01167             if (!thePdfName.Contains(theType->GetName())) continue;
01168             compStr+=thePdf->GetName();
01169             compStr+=",";
01170             compParams.add(*thePdf->getParameters(_theData));
01171           }
01172           delete compIter;
01173           if (""==compStr) continue;
01174           // do we have dataset for this cat?
01175           TString catCut=Form("%s==%s::%s",theCat->GetName(),
01176                               theCat->GetName(),theType->GetName());
01177           RooDataSet *subData=(RooDataSet*)_theData->reduce(catCut);
01178           if (subData->numEntries()<=0) {
01179             delete subData;
01180             continue;
01181           }
01182           // now plot the subData and pdf
01183           RooPlot *subFrame=theObs->frame(plotMin, plotMax, nBins);
01184           if (getControlBit("ParamsOnPlot"))
01185             doParamsOnPlot(subFrame, &compParams);
01186           RooLinkedList plotOpts;
01187           // X Error size relative to bin width, default 1.0
01188           Double_t xerrorscale = atof(readConfStr("XErrorSize", "1.0", 
01189                                                   getVarSec()));
01190           RooCmdArg xerrArg = XErrorSize(xerrorscale);
01191           plotOpts.Add((TObject*)&xerrArg);
01192           RooCmdArg datErrArg = DataError(RooAbsData::SumW2);
01193           if (subData->isWeighted()) // use sumw2
01194             plotOpts.Add((TObject*)&datErrArg);
01195           subData->plotOn(subFrame, plotOpts);
01196           _thisSimPdfWOP->plotOn(subFrame, plotCutNormCmd, Components(compStr),
01197                                  ProjWData(*subData));
01198           if (getControlBit("Chi2OnPlot")) doChi2OnPlot(subFrame);
01199           subFrame->SetNameTitle(Form("%s_%s_%s",theObs->GetName(),
01200                                       _thePdf->GetName(), catCut.Data()),
01201                                  Form("%s %s %s",theObs->GetTitle(),
01202                                       _thePdf->GetTitle(), catCut.Data()));
01203           //#ifndef RAR_USE_ROOT5
01204           subFrame->SetTitleSize(0.05);
01205           //#endif
01206           plotList.Add(subFrame);
01207           
01208           delete subData;
01209         }
01210         delete tIter;
01211       }
01212       delete cIter;
01213       delete simComps;
01214       // now plot for full simPdf
01215       {
01216         RooPlot *subFrame=theObs->frame(plotMin, plotMax, nBins);
01217         if (getControlBit("ParamsOnPlot")) doParamsOnPlot(subFrame);
01218         RooLinkedList plotOpts;
01219         Double_t xerrorscale = atof(readConfStr("XErrorSize", "1.0", 
01220                                                 getVarSec()));
01221         RooCmdArg xerrArg = XErrorSize(xerrorscale);
01222         plotOpts.Add((TObject*)&xerrArg);
01223         RooCmdArg datErrArg = DataError(RooAbsData::SumW2);
01224         if (_theData->isWeighted()) // use sumw2
01225           plotOpts.Add((TObject*)&datErrArg);
01226         _theData->plotOn(subFrame, plotOpts);
01227         _thisSimPdfWOP->plotOn(subFrame, ProjWData(*_theData),plotCutNormCmd);
01228         if (getControlBit("Chi2OnPlot")) doChi2OnPlot(subFrame);
01229         subFrame->SetNameTitle(Form("%s_%s",theObs->GetName(),
01230                                     _thisSimPdfWOP->GetName()),
01231                                Form("%s %s",theObs->GetTitle(),
01232                                     _thisSimPdfWOP->GetTitle()));
01233         //#ifndef RAR_USE_ROOT5
01234         subFrame->SetTitleSize(0.05);
01235         //#endif
01236         plotList.Add(subFrame);
01237       }
01238     }
01239     if ("simFitOnly"==pdfFitStr) continue;
01240     // now categorized plot for the total pdf
01241     TString plotWCatStr=
01242       readConfStr(Form("plotWCat_%s",theObs->GetName()),"notSet",getVarSec());
01243     if ("notSet"==plotWCatStr)
01244       plotWCatStr=readConfStr("plotWCat","no",getVarSec());
01245     if ("no"!=plotWCatStr) {
01246       rarStrParser plotCatParser=plotWCatStr;
01247       while (plotCatParser.nArgs()>0) {
01248         TString catName=plotCatParser[0];
01249         plotCatParser.Remove();
01250         RooAbsCategory *plotCat=(RooAbsCategory *)_fullObs->find(catName);
01251         if (!plotCat) {
01252           cout<<" W A R N I N G !"<<endl
01253               <<" category "<<catName<<" does not exist!"<<endl;
01254           continue;
01255         }
01256         TIterator *cIter=plotCat->typeIterator();
01257         RooCatType *cType(0);
01258         while(cType=(RooCatType*)cIter->Next()) {
01259           TString catCut=catName+"=="+catName+"::"+cType->GetName();
01260           RooDataSet *catSliceData=(RooDataSet*)_theData->reduce(catCut);
01261           RooPlot *subFrame=theObs->frame(plotMin, plotMax, nBins);
01262           if (catSliceData->isWeighted()) // use sumw2
01263             catSliceData->plotOn(subFrame, DataError(RooAbsData::SumW2));
01264           else
01265             catSliceData->plotOn(subFrame);
01266           _thePdf->plotOn(subFrame, ProjWData(*catSliceData), plotCutNormCmd);
01267           subFrame->SetNameTitle(Form("%s_%s_%s",theObs->GetName(),
01268                                       _thePdf->GetName(), catCut.Data()),
01269                                  Form("%s %s %s",theObs->GetTitle(),
01270                                       _thePdf->GetTitle(), catCut.Data()));
01271           plotList.Add(subFrame);
01272           delete catSliceData;
01273         }
01274         delete cIter;
01275       }
01276     }
01277     // pdfPlot for total pdf
01278     // get plot frame
01279     frame=theObs->frame(plotMin, plotMax, nBins);
01280     // get params
01281     if (getControlBit("ParamsOnPlot")) doParamsOnPlot(frame);
01282     // data points
01283     RooLinkedList plotOpts;
01284     Double_t xerrorscale = atof(readConfStr("XErrorSize", "1.0", getVarSec()));
01285     RooCmdArg xerrArg = XErrorSize(xerrorscale);
01286     plotOpts.Add((TObject*)&xerrArg);
01287     RooCmdArg datErrArg = DataError(RooAbsData::SumW2);
01288     if (_theData->isWeighted()) // use sumw2
01289       plotOpts.Add((TObject*)&datErrArg);
01290     _theData->plotOn(frame, plotOpts);
01291     // pdf
01292     _thePdf->plotOn(frame, ProjWData(*projWData), plotCutNormCmd);
01293     // put chi2 on
01294     if (getControlBit("Chi2OnPlot")) doChi2OnPlot(frame);
01295     // any sub plots?
01296     if (getControlBit("CompsOnPlot")) {
01297       Int_t lineStyle=4;
01298       Int_t nComp=_subPdfs.getSize();
01299       for (Int_t j=0; j<nComp; j++) {
01300         if (lineStyle>2) lineStyle=2;
01301         else lineStyle=4;
01302         RooAbsPdf *theSubPdf=(RooAbsPdf*)_subPdfs.at(j);
01303         _thePdf->plotOn(frame, ProjWData(*projWData),LineWidth(2),
01304                         plotCutNormCmd, Components(theSubPdf->GetName()),
01305                         Name(theSubPdf->GetName()),
01306                         LineStyle(lineStyle),LineColor(getColor(j)));
01307       }
01308     }
01309     // any sub data points
01310     if (getControlBit("CompsOnPlot")&&getControlBit("CompsDataOnPlot")) {
01311       RooAddPdf *refPdf(0);
01312       // do we have ref pdf in the config?
01313       TString refName=readConfStr("compsDataOnPlot", "no", getVarSec());
01314       if ("yes"==refName) refPdf=(RooAddPdf*)_thePdf;
01315       else {
01316         // ref to RooRarFit PDF?
01317         rarBasePdf *rarRef=(rarBasePdf*)(_rarPdfs.FindObject(refName));
01318         if (rarRef) refPdf=(RooAddPdf*)rarRef->getPdf();
01319         else {
01320           cout<<"Can not find ref pdf named "<<refName<<endl
01321               <<"Use the default one: "<<GetName()<<endl;
01322           refPdf=(RooAddPdf*)_thePdf;
01323         }
01324       }
01325       // make sure the refPdf is RooAddPdf
01326       if (refPdf->ClassName()!=TString("RooAddPdf")) {
01327         cout <<" The reference pdf "<<refPdf->GetName()
01328              <<" for data point plotting is not RooAddPdf"<<endl;
01329         exit(-1);
01330       }
01331       Int_t nComp=refPdf->pdfList().getSize();
01332       // now create nComp histograms for data points
01333       TList hList;
01334       for (Int_t j=0; j<nComp; j++) {
01335         TH1F *h=new
01336           TH1F(Form("subDP_%s_%s", theObs->GetName(),
01337                     refPdf->pdfList()[j].GetName()), "sub data points",
01338                frame->GetXaxis()->GetNbins(),
01339                frame->GetXaxis()->GetXmin(), frame->GetXaxis()->GetXmax());
01340         h->Sumw2();
01341         h->SetMarkerStyle(8+j);
01342         hList.Add(h);
01343       }
01344       // 
01345       refPdf->attachDataSet(*_theData);
01346       // now get data point plots
01347       Int_t nEvts=_theData->numEntries();
01348       RooArgSet *normSet=(RooArgSet *)_theData->get(0);
01349       for (Int_t evtIdx=0; evtIdx<nEvts; evtIdx++) {
01350         _theData->get(evtIdx);
01351         Double_t val=((RooAbsReal*)normSet->find(theObs->GetName()))->getVal();
01352         Double_t totProb=refPdf->getVal(normSet);
01353         Double_t lastW=1.;
01354         for (Int_t j=0; j<nComp; j++) {
01355           TH1F *h=(TH1F *)hList.At(j);
01356           if (j<nComp-1) {
01357             RooAbsPdf *subPdf=(RooAbsPdf*)refPdf->pdfList().at(j);
01358             RooAbsReal *subCoef=(RooAbsReal*)refPdf->coefList().at(j);
01359             Double_t w=subPdf->getVal(normSet)*subCoef->getVal()/totProb;
01360             lastW-=w;
01361             h->Fill(val, w);
01362           } else {
01363             h->Fill(val, lastW);
01364           }
01365         }
01366       }
01367       // add the histograms into sub frames
01368       for (Int_t j=0; j<nComp; j++) {
01369         RooPlot *subframe=theObs->frame(nBins);
01370         subframe->SetName(Form("sub_%s_%s", theObs->GetName(),
01371                                refPdf->pdfList()[j].GetName()));
01372         RooAbsPdf *theSubPdf=(RooAbsPdf*)_subPdfs.at(j);
01373         TH1F *h=(TH1F *)hList.At(j);
01374         // convert TH1 to RooHist
01375         RooHist *rHist=new RooHist(*h);
01376         subframe->addPlotable(rHist, "P");
01377         theSubPdf->plotOn(subframe, ProjWData(*projWData));
01378         plotList.Add(subframe);
01379       }
01380     }
01381     // set frame's attrs
01382     frame->SetNameTitle(Form("%s_%s", theObs->GetName(), _thePdf->GetName()),
01383                         Form("%s %s",theObs->GetTitle(),_thePdf->GetTitle()));
01384     //#ifndef RAR_USE_ROOT5
01385     frame->SetTitleSize(0.05);
01386     //#endif
01387     plotList.Add(frame);
01388   }
01389   // restore params
01390   readFromStr(*theParams, coeffParamSSaver);
01391   // restore obs saved
01392   readFromStr(_fObsSet, obsSaveStr);
01393   setControlBits("PdfPlotDone");
01394   return frame;
01395 }
01396 
01415 RooPlot *rarBasePdf::doParamsOnPlot(RooPlot* frame, RooArgSet *params,
01416                                     Int_t sigDigits, Option_t *options,
01417                                     Double_t xmin,Double_t xmax, Double_t ymax)
01418 {
01419   // parse the options
01420   TString opts = options;
01421   opts.ToLower();
01422   Bool_t showConstants= opts.Contains("c");
01423   
01424   // calculate the box's size, adjusting for constant parameters
01425   if (!params) params = _thePdf->getParameters(_theData);
01426   TIterator* pIter = params->createIterator();
01427   
01428   Int_t nLine(0);
01429   Int_t nFreeParam(0);
01430   Double_t ymin(ymax), dy(0.06);
01431   RooRealVar *var = 0;
01432   while(var=(RooRealVar*)pIter->Next()) {
01433     if (var->ClassName()!=TString("RooRealVar")) continue;
01434     if (!var->isConstant()) nFreeParam++;
01435     if(showConstants || !var->isConstant()) {
01436       ymin-= dy;
01437       nLine++;
01438     }
01439   }
01440   if(getControlBit("Chi2OnPlot")) {
01441     ymin-= dy;
01442     nLine++;
01443   }
01444   if (nLine<=0) nLine=1;
01445   Double_t step=1./nLine;
01446   
01447   // create the box and set its options
01448   TPaveText *box= new TPaveText(xmin,ymin,xmax,ymax,"BRNDC");
01449   if(!box) return frame;
01450   box->SetFillColor(0);
01451   box->SetBorderSize(1);
01452   box->SetTextAlign(12);
01453   box->SetTextSize(0.04F);
01454   box->SetFillStyle(1001);
01455   box->SetFillColor(0);
01456   TText *text(0); 
01457   pIter->Reset() ;
01458   Double_t y=1-.4*step;
01459   // add chi2 if specified
01460   if(getControlBit("Chi2OnPlot")) {
01461     text= box->AddText(0, y, Form("%d", nFreeParam));
01462     y-=step;
01463   }
01464   while(var=(RooRealVar*)pIter->Next()) {
01465     if (var->ClassName()!=TString("RooRealVar")) continue;
01466     if(var->isConstant() && !showConstants) continue;
01467     if(var->GetName()==TString(var->getPlotLabel())) {
01468       TString label=var->GetTitle();
01469       label.ReplaceAll(" ({", " _{");
01470       label.ReplaceAll(" (", " _{");
01471       label.ReplaceAll("})", "}");
01472       label.ReplaceAll(")", "}");
01473       var->setPlotLabel(label);
01474     }
01475     TString *formatted= var->format(sigDigits, opts.Data());
01476     text= box->AddText(0, y, formatted->Data());
01477     y-=step;
01478     delete formatted;
01479   }
01480   
01481   // Add box to frame 
01482   frame->addObject(box) ;
01483   
01484   delete pIter ;
01485   return frame ;
01486 }
01487 
01500 RooPlot *rarBasePdf::doChi2OnPlot(RooPlot* frame)
01501 {
01502   if (!getControlBit("ParamsOnPlot")) return frame;
01503   TPaveText* pbox = (TPaveText*) frame->findObject("TPave");
01504   if (!pbox) return frame;
01505   if (!pbox->GetSize()) return frame;
01506   TLatex *text=(TLatex*)pbox->GetLine(0);
01507   if ("nbin"==readConfStr("chi2OnPlot", "dof", getVarSec())) {
01508     text->SetTitle(Form("%s%.3f", "#chi^{2}/n = ", frame->chiSquare()));
01509   } else { // DOF
01510     Int_t nFreeParam=atoi(text->GetTitle());
01511     //text->SetLimitIndiceSize(1);
01512     //text->SetTitle(Form("%s%.3f", "#chi^{2}/_{^{DOF}} = ",
01513     text->SetTitle(Form("%s%.3f", "#chi^{2} / ndf = ",
01514                         frame->chiSquare(nFreeParam)));
01515   }
01516   
01517   return frame;
01518 }
01519 
01524 RooAbsPdf *rarBasePdf::getProtGen()
01525 {
01526   if (_theProtGen) return _theProtGen;
01527   _theProtGen=_myDummyPdf;
01528   if (_protGenPdfs.getSize()>0) {
01529     RooArgList pdfList(_protGenPdfs);
01530     pdfList.add(*_theProtGen);
01531     _theProtGen=new RooProdPdf
01532       (Form("protGen_%s",GetName()),Form("Gen for %s",GetName()), pdfList);
01533   }
01534   
01535   return _theProtGen;
01536 }
01537 
01545 RooAbsPdf *rarBasePdf::getPdfWOvar(RooArgList ignoredObs)
01546 {
01547   RooAbsPdf *thePdf=_thePdf;
01548   // check if thePdf depends on theVar
01549   if (!(thePdf->dependsOn(ignoredObs))) return thePdf;
01550   // create a dummy for it
01551   TString theName=Form("the_%s_%s", ignoredObs[0].GetName(), GetName());
01552   TString theTitle=Form("%s w/o %s", GetTitle(), ignoredObs[0].GetName());
01553   thePdf=new RooGenericPdf(theName, theTitle, "1",
01554                            *_thePdf->getParameters(_fullObs));
01555   //thePdf->Print("v");
01556   //thePdf->Print();
01557   return thePdf;
01558 }
01559 
01566 RooAbsPdf *rarBasePdf::getDPdfWvar(RooRealVar *theVar)
01567 {
01568   RooAbsPdf *retVal(0);
01569   if (!_thePdf->dependsOn(*theVar)) return retVal;
01570   retVal=_thePdf;
01571   if (_thisSimPdf) retVal=_thisSimPdf;
01572   return retVal;
01573 }
01574 
01583 void rarBasePdf::setControlBit(TString controlBitStr, TString bitConfigStr,
01584                                TString configSec)
01585 {
01586   // bit
01587   TString bit="yes";
01588   // first get bit name
01589   if (controlBitStr.BeginsWith("no")) {
01590     controlBitStr.Replace(0, 2, "");
01591     bit="no";
01592   }
01593   // then construct config name
01594   if (""==bitConfigStr) bitConfigStr=controlBitStr;
01595   if (""==configSec) configSec=getVarSec();
01596   bit=readConfStr(bitConfigStr, bit, configSec);
01597   if ("no"==bit) setControlBits("no"+controlBitStr);
01598   else setControlBits(controlBitStr);
01599 }
01600 
01611 void rarBasePdf::setControlBits(TString controlBitsStr)
01612 {
01613   rarStrParser controlBitsStrParser=controlBitsStr;
01614   for (Int_t i=0; i<controlBitsStrParser.nArgs(); i++) {
01615     TString myBit=controlBitsStrParser[i];
01616     // remove any such bit
01617     rarStrParser masterStrParser=_controlStr;
01618     _controlStr="";
01619     for (Int_t j=0; j<masterStrParser.nArgs(); j++) {
01620       if (myBit==masterStrParser[j]) continue;
01621       if ("no"+myBit==masterStrParser[j]) continue;
01622       if (myBit=="no"+masterStrParser[j]) continue;
01623       _controlStr+=" "+masterStrParser[j];
01624     }
01625     // add the bit
01626     _controlStr+=" "+myBit;
01627   }
01628 }
01629 
01638 Bool_t rarBasePdf::getControlBit(TString controlBitStr)
01639 {
01640   rarStrParser controlStrParser=_controlStr;
01641   for (Int_t i=0; i<controlStrParser.nArgs(); i++) {
01642     if (controlBitStr==controlStrParser[i]) return kTRUE;
01643     if ("no"+controlBitStr==controlStrParser[i]) return kFALSE;
01644     if (controlBitStr=="no"+controlStrParser[i]) return kFALSE;
01645   }
01646   // not exist
01647   if (controlBitStr.BeginsWith("no")) return kTRUE;
01648   return kFALSE;
01649 }
01650 
01654 Int_t rarBasePdf::getColor(Int_t i)
01655 {
01656   return _rarColors[i%NCOLORS];
01657 }

Generated on 30 Oct 2013 for RooRarFit by  doxygen 1.4.7