RooRarFit: rarMLPdf.cc Source File

00001 /*****************************************************************************
00002 * Project: BaBar detector at the SLAC PEP-II B-factory
00003 * Package: RooRarFit
00004  *    File: $Id: rarMLPdf.cc,v 1.27 2011/06/16 13:18:50 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 ML model class for RooRarFit
00014 //
00015 // BEGIN_HTML
00016 // This class provides ML model class for RooRarFit
00017 // END_HTML
00018 //
00019 
00020 #include "RooRarFit/rarVersion.hh"
00021 
00022 #include "Riostream.h"
00023 
00024 #include "RooFitCore/Roo1DTable.hh"
00025 #include "RooFitCore/RooAbsPdf.hh"
00026 #include "RooFitCore/RooArgList.hh"
00027 #include "RooFitCore/RooDataSet.hh"
00028 #include "RooFitCore/RooFitResult.hh"
00029 #include "RooFitCore/RooMCStudy.hh"
00030 #include "RooFitCore/RooProdPdf.hh"
00031 #include "RooFitCore/RooRealVar.hh"
00032 #include "RooFitCore/RooStringVar.hh"
00033 
00034 #include "RooFitCore/RooAddPdf.hh"
00035 
00036 #include "RooRarFit/rarMLFitter.hh"
00037 #include "RooRarFit/rarMLPdf.hh"
00038 
00039 ClassImp(rarMLPdf)
00040   ;
00041 
00045 rarMLPdf::rarMLPdf()
00046   : rarAdd(),
00047     _specialStr("")
00048 {
00049   init();
00050 }
00051 
00064 rarMLPdf::rarMLPdf(const char *configFile, const char *configSec,
00065                    const char *configStr,
00066                    rarDatasets *theDatasets, RooDataSet *theData,
00067                    const char *name, const char *title)
00068   : rarAdd(configFile, configSec, configStr,
00069               theDatasets, theData, name, title, kFALSE, kFALSE),
00070     _specialStr("")
00071 {
00072   init();
00073 }
00074 
00075 rarMLPdf::~rarMLPdf()
00076 {
00077 }
00078 
00104 void rarMLPdf::init()
00105 {
00106   // make sure it is extended
00107   if (_nCoeff!=_nComp) {
00108     cout<<"Pdf model should be extended pdf"<<endl;
00109     exit(-1);
00110   }
00111   // add compCat type, the last type for embedded events
00112   for (Int_t i=0; i<=_nComp; i++) _compCat.defineType(Form("%02d",i));
00113   // save the original coeff List anyway
00114   _sCoeffs.add(_coeffs);
00115   // to deal with simu fit and splitting of coeffs
00116   TString masterSec=getMasterSec();
00117   Bool_t simFit=getControlBit("SimFit")&&
00118     ("simple"!=readConfStr("simultaneousFit", "yes",masterSec));
00119   // for split specials
00120   if (simFit) {
00121     createAbsVars("splitSpecials", &_specialSet);
00122     // run master sec's
00123     TString myVarSec=getVarSec();
00124     setVarSec(masterSec);
00125     createAbsVars("splitSpecials", &_specialSet);
00126     setVarSec(myVarSec);
00127   }
00128   TString yieldSplitMethod=readConfStr("yieldSplitMethod", "auto", masterSec);
00129   if (simFit && ("manual"!=yieldSplitMethod)) {
00130     // use master sec config
00131     TString myVarSec=getVarSec();
00132     setVarSec(masterSec);
00133     // get coeff list, coeff param set, param list
00134     RooArgSet coeffParamSet("Coeff Parameter Set");
00135     for (Int_t i=0; i<_nCoeff; i++) {
00136       coeffParamSet.add(*(_coeffs[i].getParameters(*_fullObs)));
00137     }
00138     RooArgList coeffParamList(coeffParamSet);
00139     // let's make sure coeffs are not split in the mlFitter
00140     // coz we are going to split them here
00141     TString pdfSplitStr=readConfStr(Form("the_%s", GetName()), "", masterSec);
00142     if ("semi"==yieldSplitMethod) {
00143       for (Int_t i=0; i<_coeffs.getSize(); i++) {
00144         if (pdfSplitStr.Contains(_coeffs[i].GetName())) {
00145           cout<<"Please do not include "<<_coeffs[i].GetName()
00146               <<" or its dependents in your simPdfBuilder configs"<<endl
00147               <<"in section \""<<masterSec<<"\":"<<endl
00148               <<pdfSplitStr<<endl
00149               <<"We will split it for you."<<endl;
00150           exit(-1);
00151         }
00152       }
00153       for (Int_t i=0; i<coeffParamList.getSize(); i++) {
00154         if (pdfSplitStr.Contains(coeffParamList[i].GetName())) {
00155           cout<<"Please do not include "<<coeffParamList[i].GetName()
00156               <<" or its dependents in your simPdfBuilder configs"<<endl
00157               <<"in section \""<<masterSec<<"\":"<<endl
00158               <<pdfSplitStr<<endl
00159               <<"We will split it for you."<<endl;
00160           exit(-1);
00161         }
00162       }
00163     }
00164     // get splitting cats
00165     TString physCat=getFitter()->getPhysCat();
00166     TString splitCats=getFitter()->getSplitCats();
00167     // build split coeffs
00168     _coeffs.removeAll();
00169     for (Int_t i=0; i<_nCoeff; i++) {
00170       // formula for the final coeff
00171       TString coeffArgsStr=_sCoeffs[i].GetName();
00172       TString formulaStr="@0";
00173       // get frac split rules
00174       TString fracRuleStr=readConfStr("fracRule", splitCats, masterSec);
00175       fracRuleStr=readConfStr("fracRule_"+coeffArgsStr,fracRuleStr, masterSec);
00176       rarStrParser catsStrParser=fracRuleStr;
00177       Int_t nSCat=catsStrParser.nArgs();
00178       // create multiple cat coeffs
00179       RooArgSet splitCatSet(getFitter()->getSplitCatSet());
00180       for (Int_t j=0; j<nSCat; j++) {
00181         TString catName=catsStrParser[j];
00182         RooAbsCategory *theCat=getFitter()->getSplitCat(splitCatSet, catName);
00183         if (!theCat) {
00184           cout<<"Can not find cat "<<catName<<endl
00185               <<"Please make sure the cat you are using exists"<<endl;
00186           exit(-1);
00187         }
00188         // rename the catName
00189         catName=theCat->GetName();
00190         // cat ceoff name
00191         TString catCoeffName=Form("Frac_%s_%s",
00192                                   _sCoeffs[i].GetName(), catName.Data());
00193         TString catCoeffStr=readConfStr(catCoeffName, "notSet", masterSec);
00194         if ("notSet"==catCoeffStr)
00195           catCoeffStr=Form("%s %s RooRealVar \"cat coeff\" C 1. 0. 1.",
00196                            catCoeffName.Data(), catCoeffName.Data());
00197         else
00198           catCoeffStr=catCoeffName+" "+catCoeffStr;
00199         RooRealVar *catCoeff=(RooRealVar*)createAbsVar(catCoeffStr);
00200         catCoeffName=catCoeff->GetName();
00201         // save the var name
00202         saveFracName(catCoeffName);
00203         // addToConfStr("Ignored", catCoeff->GetName(), myVarSec);
00204         coeffArgsStr+=Form(" %s ", catCoeff->GetName());
00205         formulaStr+=Form("*@%d",j+1);
00206         // let's make sure catCoeff is created as RooRealVar
00207         if ("RooRealVar"!=TString(catCoeff->ClassName())) {
00208           cout<<catCoeffName<<" is created in config file as "
00209               <<catCoeff->ClassName()<<","<<endl
00210               <<"but it should be RooRealVar"<<endl;
00211           exit(-1);
00212         }
00213         // now, we create split catCoeffs for this catCoeff
00214         // first we find out how many types for this cat
00215         TString catTypesStr="";
00216         TIterator* catTypeIter = theCat->typeIterator();
00217         RooCatType *theType(0);
00218         while(theType=(RooCatType*)catTypeIter->Next()) {
00219           catTypesStr+=Form(" %s ", theType->GetName());
00220         }
00221         delete catTypeIter;
00222         //cout<<catTypesStr<<endl;
00223         rarStrParser catTypesStrParser=catTypesStr;
00224         if (catTypesStrParser.nArgs()<1) continue; // nothing in the cat
00225         TString cat1FormulaStr="1.";
00226         TString cat1ArgSet="";
00227         Roo1DTable* theTable(0);
00228         if (_theData) theTable = _theData->table(*theCat);
00229         // check if we have fracSrc config
00230         Bool_t useFracSrc=kFALSE;
00231         TString fracSrcStr=
00232           readConfStrCnA(Form("fracSrc_%s_%s", _sCoeffs[i].GetName(),
00233                               catName.Data()),"notSet");
00234         if ("notSet"==fracSrcStr)
00235           fracSrcStr=
00236             readConfStrCnA(Form("fracSrc_%s",_sCoeffs[i].GetName()), "notSet");
00237         if ("notSet"==fracSrcStr)
00238           fracSrcStr=readConfStrCnA("fracSrc", "notSet");
00239         if ("notSet"!=fracSrcStr) {
00240           RooDataSet *theData=getDatasets()->getData(fracSrcStr);
00241           if (theData) {
00242             theTable = theData->table(*theCat);
00243             useFracSrc=kTRUE;
00244             cout<<" Using table from "<<theData->GetName()
00245                 <<" for "<<_sCoeffs[i].GetName()<<" "<<catName<<endl;
00246             theTable->Print();
00247           }
00248         }
00249         // now we need to find out which type to be the non-free
00250         Int_t nonIdx=0;
00251         // first for auto mode
00252         if ("auto"==yieldSplitMethod) {
00253           Bool_t notFound=(catTypesStrParser.nArgs()>1); // if #types>1
00254           for (Int_t k=0; k<catTypesStrParser.nArgs(); k++) {
00255             TString myNonFreeName=catCoeffName+"["+catTypesStrParser[k]+"]";
00256             Int_t strIdx=pdfSplitStr.Index(myNonFreeName);
00257             if (strIdx<0) continue;
00258             // check cat seq
00259             TString subpdfSplitStr=pdfSplitStr(0,strIdx);
00260             Int_t colonIdx=subpdfSplitStr.Last(':');
00261             if (colonIdx<=0) break;
00262             rarStrParser pdfSplitStrParser=TString(subpdfSplitStr(0,colonIdx));
00263             //cout<<pdfSplitStr<<endl<<myNonFreeName<<endl<<strIdx<<endl;
00264             // get cats
00265             TString splitCatName=
00266               pdfSplitStrParser[pdfSplitStrParser.nArgs()-1];
00267             splitCatName.ReplaceAll(",", "_");
00268             //cout<<splitCatName<<endl;
00269             if (splitCatName==catName) {
00270               nonIdx=k;
00271               notFound=kFALSE;
00272             }
00273             break;
00274           }
00275           if (notFound) {
00276             cout<<" Can not find splitting specif. for"<<endl
00277                 <<" "<<catCoeffName<<" wrt types:"<<endl
00278                 <<" "<<catTypesStr<<endl
00279                 <<" in the_"<<GetName()<<endl<<pdfSplitStr<<endl
00280                 <<" from section ["<<getVarSec()<<"]"<<endl;
00281             exit(-1);
00282           }
00283         } else { // semi mode
00284           TString nonFreeStr=readConfStr("nonFreeCatTypes", "", getVarSec());
00285           for (Int_t k=0; k<catTypesStrParser.nArgs(); k++) {
00286             if (nonFreeStr.
00287                 Contains(catCoeffName+"["+catTypesStrParser[k]+"]")) {
00288               nonIdx=k;
00289               break;
00290             }
00291           }
00292         }
00293         for (Int_t k=0; k<catTypesStrParser.nArgs(); k++) {
00294           // make sure the name is a cat type
00295           const RooCatType *theType=theCat->lookupType(catTypesStrParser[k]);
00296           if (!theType) {
00297             cout <<"The cat type "<<catTypesStrParser[k]
00298                  <<" does not exist in cat "<<theCat->GetName()<<endl;
00299             exit(-1);
00300           }
00301           if (k==nonIdx) continue;
00302           TString catCoeffTypeName=catCoeffName+"_"+catTypesStrParser[k]+"";
00303           Double_t typeFrac=1./catTypesStrParser.nArgs();
00304           if (theTable) typeFrac=theTable->getFrac(catTypesStrParser[k]);
00305           RooAbsReal *catCoeffType=(RooAbsReal*)
00306             createAbsVar(Form("%s %s RooRealVar \"cat type coeff\" C %f 0 1",
00307                               catCoeffTypeName.Data(), catCoeffTypeName.Data(),
00308                               typeFrac));
00309           // save the var name
00310           saveFracName(catCoeffTypeName);
00311           if (useFracSrc)
00312             addToConfStr("Ignored", catCoeffType->GetName(), myVarSec);
00313           _specialSet.add(*catCoeffType);
00314           if (2==catTypesStrParser.nArgs()) _asymSet.add(*catCoeffType);
00315           Int_t fvIdx = k<nonIdx ? k : k-1;
00316           cat1FormulaStr+=Form("-@%d", fvIdx);
00317           cat1ArgSet+=Form(" %s ", catCoeffType->GetName());
00318         }
00319         // for nonfree cat type
00320         TString catCoeffTypeName=catCoeffName+"_"+catTypesStrParser[nonIdx]+"";
00321         // old way, ie RoorarFit way, RooRealVar
00322         RooAbsReal *catCoeffType=(RooAbsReal*)
00323           createAbsVar(Form("%s %s RooFormulaVar %s %s",
00324                             catCoeffTypeName.Data(), catCoeffTypeName.Data(),
00325                             cat1FormulaStr.Data(), cat1ArgSet.Data()));
00326         // make sure it is RooFormulaVar
00327         if ("RooFormulaVar"!=TString(catCoeffType->ClassName())) {
00328           cout<<catCoeffTypeName<<" is created in config file as "
00329               <<catCoeffType->ClassName()<<","<<endl
00330               <<"but it should be RooFormulaVar"<<endl;
00331           exit(-1);
00332         }
00333         _specialSet.add(*catCoeffType);
00334         if (catTypesStrParser.nArgs()>1) _cat1Set.add(*catCoeffType);
00335         if ("semi"==yieldSplitMethod) {
00336           _specialStr+=Form(" %s : %s ", catName.Data(), catCoeffName.Data());
00337         }
00338       }
00339       // do we still have any left over?
00340       if (splitCatSet.getSize()>0) {
00341         cout<<" No splitting rules for "<<_sCoeffs[i].GetName()<<" wrt:"<<endl;
00342         splitCatSet.Print();
00343         cout<<" Please check config fracRule or fracRule_"
00344             <<_sCoeffs[i].GetName()<<" in section \""<<masterSec<<"\":"<<endl;
00345         exit(-1);
00346       }
00347       TString coeffName=Form("theSim_%s", _sCoeffs[i].GetName());
00348       //cout<<coeffName<<endl
00349       //<<formulaStr<<endl
00350       //<<coeffArgsStr<<endl;
00351       RooAbsReal *coeff=(RooAbsReal*)
00352         createAbsVar(Form("%s %s RooFormulaVar %s %s",
00353                           coeffName.Data(), coeffName.Data(),
00354                           formulaStr.Data(), coeffArgsStr.Data()));
00355       //coeff->Print("v");
00356       _coeffs.add(*coeff);
00357     }
00358     //_specialSet.Print("v");
00359     // restore myVarSec
00360     setVarSec(myVarSec);
00361   }
00362   // build pdf function
00363   _thePdf=new RooAddPdf(Form("the_%s", GetName()), _pdfType+" "+GetTitle(),
00364                         _subPdfs, _coeffs);
00365   //_thePdf->Print("v");
00366   cout<<"done init of rarMLPdf for "<<GetName()<<endl<<endl;
00367 }
00368 
00375 RooArgSet rarMLPdf::getSpecialSet(TString setName)
00376 {
00377   if ("cat1Set"==setName) return _cat1Set;
00378   if ("asymSet"==setName) return _asymSet;
00379   return _specialSet; // default
00380 }
00381 
00386 RooAbsPdf *rarMLPdf::getProtGen()
00387 {
00388   if (_theProtGen) return _theProtGen;
00389   // call super class' getProtGen
00390   rarCompBase::getProtGen();
00391   // construct prototype var generator from components
00392   RooArgList compList;
00393   for (Int_t i=0; i<_nComp; i++) {
00394     if (_protGenPdfs.getSize()>0) {
00395       RooArgList pdfList(_protGenPdfs);
00396       pdfList.add(_subProtGenPdfs[i]);
00397       TString compName=
00398         Form("prodProtGen_%s_%s",GetName(),_subPdfs[i].GetName());
00399       TString compTitle=
00400         Form("Prod Prot Gen %s %s",GetName(),_subPdfs[i].GetName());
00401       RooAbsPdf *thePdf=new RooProdPdf (compName, compTitle, pdfList);
00402       compList.add(*thePdf);
00403     } else compList.add(_subProtGenPdfs[i]);
00404   }
00405   _theProtGen=new RooAddPdf
00406     (Form("gen_%s", GetName()), Form("Gen Add for %s", GetName()),
00407      compList, _coeffs);
00408   
00409   cout<<"getProtGen: Expected number of events "
00410       <<_theProtGen->expectedEvents(0)<<endl;
00411   
00412   cout<<"protGen created for "<<GetName()<<endl;
00413   _theProtGen->Print();
00414   //_theProtGen->Print("v");
00415   return _theProtGen;
00416 }
00417 
00428 RooPlot *rarMLPdf::doPdfPlot(TList &plotList, TString pdfList)
00429 {
00430   return rarCompBase::doPdfPlot(plotList, pdfList);
00431 }