RooRarFit: rarConfig.cc Source File

00001 /*****************************************************************************
00002 * Project: BaBar detector at the SLAC PEP-II B-factory
00003 * Package: RooRarFit
00004  *    File: $Id: rarConfig.cc,v 1.61 2012/12/07 11:46:38 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 config class for RooRarFit
00014 //
00015 // BEGIN_HTML
00016 // This class provides config class for RooRarFit
00017 // END_HTML
00018 //
00019 
00020 #include "RooRarFit/rarVersion.hh"
00021 
00022 #include "Riostream.h"
00023 #include <sstream>
00024 #include <vector>
00025 
00026 #include <ctype.h>
00027 
00028 #include "TFile.h"
00029 #include "TTree.h"
00030 #include "TSystem.h"
00031 #include "TMatrixDSym.h"
00032 
00033 #include "RooFitCore/Roo1DTable.hh"
00034 #include "RooFitCore/RooAbsCategory.hh"
00035 #include "RooFitCore/RooArgList.hh"
00036 #include "RooFitCore/RooCategory.hh"
00037 #include "RooFitCore/RooConstVar.hh"
00038 #include "RooFitCore/RooDataSet.hh"
00039 #include "RooFitCore/RooDataHist.hh"
00040 #include "RooFitCore/RooMappedCategory.hh"
00041 #include "RooFitCore/RooRandom.hh"
00042 #include "RooFitCore/RooRealVar.hh"
00043 #include "RooFitCore/RooStringVar.hh"
00044 #include "RooFitCore/RooThresholdCategory.hh"
00045 #include "RooFitCore/RooGlobalFunc.hh"
00046 
00047 #include "RooFitCore/RooSuperCategory.hh"
00048 #include "RooFitModels/RooUnblindOffset.hh"
00049 #include "RooFitModels/RooUnblindPrecision.hh"
00050 
00051 #include "RooRarFit/rarAdd.hh"
00052 #include "RooRarFit/rarArgusBG.hh"
00053 #include "RooRarFit/rarDecay.hh"
00054 #include "RooRarFit/rarBasePdf.hh"
00055 #include "RooRarFit/rarBifurGauss.hh"
00056 #include "RooRarFit/rarCBShape.hh"
00057 #include "RooRarFit/rarDatasets.hh"
00058 #include "RooRarFit/rarExp.hh"
00059 #include "RooRarFit/rarGaussModel.hh"
00060 #include "RooRarFit/rarGaussian.hh"
00061 #include "RooRarFit/rarOsipDisc.hh"
00062 #include "RooRarFit/rarGeneric.hh"
00063 #include "RooRarFit/rarHistPdf.hh"
00064 #include "RooRarFit/rarKeys.hh"
00065 #include "RooRarFit/rarMLPdf.hh"
00066 #include "RooRarFit/rarNovosibirsk.hh"
00067 #include "RooRarFit/rarBallack.hh"
00068 #include "RooRarFit/rarPoly.hh"
00069 #include "RooRarFit/rarProd.hh"
00070 #include "RooRarFit/rarSimPdf.hh"
00071 #include "RooRarFit/rarStep.hh"
00072 #include "RooRarFit/rarBinned.hh"
00073 #include "RooRarFit/rarTriGauss.hh"
00074 #include "RooRarFit/rarTwoGauss.hh"
00075 #include "RooRarFit/rarUsrPdf.hh"
00076 #include "RooRarFit/rarCruijff.hh"
00077 #include "RooRarFit/rarMultPdf.hh"
00078 #include "RooRarFit/rarLass.hh"
00079 #include "RooRarFit/rarRelBreitWigner.hh"
00080 #include "RooRarFit/rarVoigtian.hh"
00081 #include "RooRarFit/rarGounarisSakurai.hh"
00082 #include "RooRarFit/rarFlatte.hh"
00083 #include "RooRarFit/rarUniform.hh"
00084 #include "RooRarFit/rarThreshold.hh"
00085 
00086 #include "RooRarFit/rarConfig.hh"
00087 
00088 ClassImp(rarConfig)
00089   ;
00090 
00091 TList rarConfig::_rarPdfs;
00092 TList rarConfig::_rarVars;
00093 TList rarConfig::_rarOVars;
00094 RooArgSet rarConfig::_rarCats;
00095 TString rarConfig::_masterSec;
00096 TString rarConfig::_runSec;
00097 
00101 rarConfig::rarConfig()
00102   : TNamed ("",""),
00103     _configFile(""),_configSec(""), _configStr(""),
00104     _createFundamental(kFALSE), _fullObs(0)
00105 {
00106   init();
00107 }
00108 
00119 rarConfig::rarConfig(const char *configFile, const char *configSec,
00120                          const char *configStr,
00121                          const char *name, const char *title)
00122   : TNamed(name, title),
00123     _configFile(configFile), _configSec(configSec), _configStr(configStr),
00124     _createFundamental(kFALSE), _fullObs(0)
00125 {
00126   init();
00127 }
00128 
00130 rarConfig::~rarConfig()
00131 {
00132 }
00133 
00140 void rarConfig::init()
00141 {
00142   cout<<endl
00143       <<GetName()<<":"<<endl
00144       <<" init of rarConfig"<<endl
00145       <<" configSec: "<<_configSec<<endl
00146       <<" configStr: "<<_configStr<<endl;
00147   
00148   // default name schema
00149   _fullNameSchema="prefix";
00150   
00151   // check if the config file is readable
00152   ifstream ifs(_configFile);
00153   if(ifs.fail()) {
00154     cout<<"rarConfig::init() can not open config file "<<_configFile<<endl;
00155     exit(-1);
00156   }
00157   ifs.close();
00158 }
00159 
00172 TString rarConfig::getFullFileName(const TString dir, const TString aType, const TString name,
00173                                    const TString dsName, const TString msName,
00174                                    const TString cfName)
00175 {
00176   TString fileName=dir;
00177   // check dir
00178   if (gSystem->AccessPathName(dir)&&gSystem->mkdir(dir)) {
00179     cout << "getFullFileName: Can not access dir "<<dir<<endl;
00180     exit(-1);
00181   }
00182   // add configFile
00183   TString configFile=cfName;
00184   if (""==configFile) configFile=_configFile;
00185   configFile=gSystem->BaseName(configFile);
00186   // remove .config
00187   if (configFile.EndsWith(".config"))
00188     configFile.Remove(configFile.Length()-7, 7);
00189   fileName+="/"+getFileName(configFile);
00190 
00191   // add dataset sec
00192   TString datasetName = dsName;
00193   if (""==datasetName) {
00194     if (getDatasets()) datasetName=getDatasets()->getVarSec();
00195     else datasetName=getVarSec();
00196   }
00197   fileName+=getFileName(datasetName);
00198 
00199   // add master sec
00200   TString masterName = msName;
00201   if (""==masterName) masterName=getMasterSec();
00202   fileName+=getFileName(masterName);
00203 
00204   // add action type
00205   fileName+=getFileName(aType);
00206 
00207   // add action name
00208   fileName+=getFileName(name);
00209   // remove the last dot
00210   if (fileName.EndsWith(".")) fileName.Remove(fileName.Last('.'), 1);
00211   
00212   return fileName;
00213 }
00214 
00222 TString rarConfig::getFileName(const TString in_name) const
00223 {
00224   TString name = in_name;
00225   if (("none"==name)||("no"==name)||("yes"==name)||("default"==name))
00226     return "";
00227   
00228   for (Int_t i=0; i<name.Length(); i++)
00229     if (!isalnum(name[i])) name[i]='_';
00230   name+=".";
00231   
00232   return name;
00233 }
00234 
00246 TString rarConfig::readConfStr(const char *name, const char *val,
00247                                const char *secName)
00248 {
00249   // construct section name
00250   TString configSec=_configSec;
00251   if (secName) configSec=secName;
00252   TString secVarName=Form("%s_%s", name, configSec.Data());
00253   // first check if it has been read in
00254   RooStringVar *theStr=(RooStringVar*)(_configStrSet.find(secVarName));
00255   if (theStr) return theStr->getVal();
00256   // not read in yet, try to get it from config file
00257   RooArgSet strList("Read Config String List");
00258   RooStringVar strVar(name, "config string", val, 40960);
00259   strList.add(strVar);
00260   strList.readFromFile(_configFile, 0, configSec);
00261   
00262   Int_t strLen=strlen(strVar.getVal())+1;
00263   if (strLen<1024) strLen=1024;
00264   if (strVar.getVal()!=TString(val)) _configStrSet.addOwned
00265     (*(new RooStringVar(secVarName, secVarName, strVar.getVal(), strLen)));
00266   
00267   return strVar;
00268 }
00269 
00276 void rarConfig::setConfStr(const char*name, const char*val, const char*secName)
00277 {
00278   // construct section name
00279   TString configSec=_configSec;
00280   if (secName) configSec=secName;
00281   TString secVarName=Form("%s_%s", name, configSec.Data());
00282   // first check if it has been read in
00283   RooStringVar *theStr=(RooStringVar*)(_configStrSet.find(secVarName));
00284   if (!theStr) { // create it
00285     Int_t strLen=strlen(val)+1;
00286     if (strLen<1024) strLen=1024;
00287     theStr=new RooStringVar(secVarName, secVarName, val, strLen);
00288     _configStrSet.addOwned(*theStr);
00289   }
00290   theStr->setVal(val);
00291   if (!val) { // we should remove this config
00292     _configStrSet.remove(*theStr);
00293     delete theStr;
00294   }
00295 }
00296 
00303 void rarConfig::addToConfStr(const char*name,const char*val,const char*secName)
00304 {
00305   TString configStr=readConfStr(name, "", secName);
00306   setConfStr(name, configStr+" "+val, secName);
00307 }
00308 
00316 TString rarConfig::readConfStrCnA(TString configStr, TString defVal)
00317 {
00318   TString retVal=readConfStr(configStr, "notSet", getVarSec());
00319   if ("notSet"==retVal) retVal=readConfStr(configStr, defVal, _runSec);
00320   return retVal;
00321 }
00322 
00331 Bool_t rarConfig::isNumber(TString numStr)
00332 {
00333   Bool_t retVal(kTRUE);
00334   char c=numStr[0];
00335   if (isdigit(c)) return retVal;
00336   if (0!=atof(numStr)) return retVal;
00337   if ('.'==c) return retVal;
00338   if (('+'==c)||('-'==c)) return retVal;
00339   
00340   retVal=kFALSE;
00341   return retVal;
00342 }
00343 
00351 Bool_t rarConfig::isVarType(TString typeStr)
00352 {
00353   if (
00354       "RooRealVar"==typeStr ||
00355       "RooConstVar"==typeStr ||
00356       "RooUnblindOffset"==typeStr ||
00357       "RooUnblindPrecision"==typeStr ||
00358       "RooCategory"==typeStr ||
00359       "RooMappedCategory"==typeStr ||
00360       "RooThresholdCategory"==typeStr ||
00361       "RooSuperCategory"==typeStr ||
00362       "RooStringVar"==typeStr ||
00363       "RooFormulaVar"==typeStr
00364       ) return kTRUE;
00365   return kFALSE;
00366 }
00367 
00373 void rarConfig::writeToStr(RooArgSet &aSet, string &aStr)
00374 {
00375   stringstream aSaver;
00376   aSet.writeToStream((ostream&)aSaver, kFALSE);
00377   aStr=aSaver.str();  
00378 }
00379 
00385 void rarConfig::readFromStr(RooArgSet &aSet, string &aStr)
00386 {
00387   stringstream aSaver;
00388   aSaver.str(aStr);
00389   aSet.readFromStream((istream&)aSaver, kFALSE);
00390 }
00391 
00406 rarStrParser rarConfig::getVarTNTU(TString varStr, TString option,
00407                                    TString *fullName, TString *fullTitle,
00408                                    TString *varType, TString *Name,
00409                                    TString *Title, TString *Unit)
00410 {
00411   // get the parsers
00412   rarStrParser varStrParser=varStr;
00413   rarStrParser optionParser=option;
00414   TString myFullName(""), myFullTitle(""), myVarType(""),
00415     myName(""), myTitle(""), myUnit("");
00416   Bool_t fNameSpecified(kFALSE);
00417   // first for short name
00418   if (varStrParser.nArgs()>0) {
00419     myName=varStrParser[0];
00420     varStrParser.Remove();
00421   }
00422   if (varStrParser.nArgs()>0) {
00423     // check if the next token is name or not
00424     if (!isVarType(varStrParser[0])) { // it is full name, or new schema
00425       if ((varStrParser.nArgs()>1)&&isVarType(varStrParser[1])){// the old way
00426         myFullName=varStrParser[0];
00427         varStrParser.Remove();
00428         fNameSpecified=kTRUE; // full name specified
00429       } else if (varStrParser.nArgs()==1) { // full name only
00430         myFullName=varStrParser[0];
00431         varStrParser.Remove();
00432         fNameSpecified=kTRUE; // full name specified
00433       } else { // the new schema
00434         myVarType="nsRooRealVar"; // new schema RRV
00435         myFullName=myName;
00436         if (optionParser.Have("N")) {
00437           myFullName=optionParser[optionParser.Index("N")+1];
00438           fNameSpecified=kTRUE; // full name specified
00439         }
00440         if (optionParser.Have("T"))
00441           myTitle=optionParser[optionParser.Index("T")+1];
00442         if (optionParser.Have("U"))
00443           myUnit=optionParser[optionParser.Index("U")+1];
00444         // remember current index 0 for N, 1 for T and 2 for U
00445         Int_t cvIdx(0);
00446         while((varStrParser.nArgs()>0)&&(!isNumber(varStrParser[0]))) {
00447           TString flag=varStrParser[0];
00448           varStrParser.Remove();
00449           if (varStrParser.nArgs()>0) {
00450             if ("N"==flag) {
00451               myFullName=varStrParser[0];
00452               varStrParser.Remove();
00453               fNameSpecified=kTRUE; // full name specified
00454             } else if ("T"==flag) {
00455               myTitle=varStrParser[0];
00456               varStrParser.Remove();
00457             } else if ("U"==flag) {
00458               myUnit=varStrParser[0];
00459               varStrParser.Remove();
00460             } else { // OK flag is the one
00461               if (0==cvIdx) {
00462                 myFullName=flag;
00463                 fNameSpecified=kTRUE; // full name specified
00464               } else if(1==cvIdx) myTitle=flag;
00465               else if(2==cvIdx) myUnit=flag;
00466             }
00467           } else {
00468             if (0==cvIdx) {
00469               myFullName=flag;
00470               fNameSpecified=kTRUE; // full name specified
00471             } else if(1==cvIdx) myTitle=flag;
00472             else if(2==cvIdx) myUnit=flag;
00473           }
00474           cvIdx++;
00475         }
00476       }
00477     }
00478     // let's make sure the name is not notSet
00479     if ("notSet"==myName) {
00480       cout<<"var name should not be notSet, please give the var a name"<<endl
00481           <<"varStr: "<<varStr<<endl;
00482       exit(-1);
00483     }
00484     // now, find varType
00485     if ((varStrParser.nArgs()>0)&&(""==myVarType)) {
00486       if (!isVarType(varStrParser[0])) { // wrong
00487         cout<<"This is supposed to be a varType "<<varStrParser[0]<<endl;
00488         exit(-1);
00489       }
00490       myVarType=varStrParser[0];
00491       varStrParser.Remove();
00492       // then the next token is title
00493       myTitle=varStrParser[0];
00494       varStrParser.Remove();
00495     }
00496   }
00497   if (!fNameSpecified) {
00498     myFullName=getFullVarName(myName);
00499     if (myFullName==myName) myFullTitle=myTitle;
00500     else myFullTitle=myTitle; //+" "+GetTitle();
00501   } else {
00502     myFullTitle=myTitle;
00503   }
00504   
00505   if (fullName) *fullName=myFullName;
00506   if (fullTitle) *fullTitle=myFullTitle;
00507   if (varType) *varType=myVarType;
00508   if (Name) *Name=myName;
00509   if (Title) *Title=myTitle;
00510   if (Unit) *Unit=myUnit;
00511   return varStrParser;
00512 }
00513 
00521 TString rarConfig::getFullVarName(TString nameStr)
00522 {
00523   TString myName=GetName();
00524   TString myFullName;
00525   if ("prefix"==_fullNameSchema) myFullName=myName+"_"+nameStr;
00526   else if ("suffix"==_fullNameSchema) myFullName=nameStr+"_"+myName;
00527   else myFullName=nameStr;
00528   return myFullName;
00529 }
00530 
00540 TString rarConfig::getFullVarName(TString nameStr,
00541                                   TString varStr, TString option)
00542 {
00543   TString myFullName;
00544   if ((""==varStr)&&(""==option)) return getFullVarName(nameStr);
00545   getVarTNTU(varStr, option, &myFullName);
00546   return myFullName;
00547 }
00548 
00588 RooAbsReal *rarConfig::createAbsReal(const char *name, const char *title,
00589                                      const Double_t val,
00590                                      const Double_t min, const Double_t max,
00591                                      const char *unit)
00592 {
00593   RooAbsReal *theVar(0);
00594   if (max<min) {
00595     cout<<"max="<<max<<" must >= "<<"min="<<min<<endl;
00596     exit(-1);
00597   }
00598   // check option to see if it is for RooConstVar
00599   TString option=unit;
00600   Bool_t forRooConstVar(kFALSE);
00601   if (option.Contains("RooConstVar")) {
00602     forRooConstVar=kTRUE;
00603     option.ReplaceAll("RooConstVar", "");
00604   }
00605   
00606   // create the construct string
00607   TString varStr=Form("%s RooRealVar \"%s\"", name, title);
00608   if (forRooConstVar) varStr=Form("%s RooConstVar \"%s\"", name, title);
00609   if ((0==min) && (0==max)) {
00610     varStr=Form("%s %f", varStr.Data(), val);
00611   } else {
00612     varStr=Form("%s %f %f %f", varStr.Data(), val, min, max);
00613   }
00614   if (""!=option) varStr=Form("%s \"%s\"", varStr.Data(), option.Data());
00615   TString varConfigString=readConfStr(name, "notSet", getVarSec());
00616   if ("notSet"!=varConfigString) {
00617     if (isNumber(varConfigString)&&forRooConstVar) {
00618       varStr=
00619         Form("%s RooConstVar \"%s\" %s", name, title, varConfigString.Data());
00620     } else {
00621       // construct new varStr
00622       varStr=Form("%s %s", name, varConfigString.Data());
00623     }
00624   }
00625   // whatever in option right now is actually unit
00626   if (""!=option) option ="U "+option;
00627   // now add title into option
00628   option.Append(Form(" T %s ", title));
00629   // the actual function to create variable
00630   theVar=(RooAbsReal*)createAbsVar(varStr, option);
00631   return theVar;
00632 }
00633 
00648 RooAbsReal *rarConfig::createAbsReal(const char *name, const char *title,
00649                                        const Double_t val, const char *unit)
00650 {
00651   return createAbsReal(name, title, val, 0, 0, unit);
00652 }
00653 
00708 RooAbsArg *rarConfig::createAbsVar(const char *varStr, const char *option)
00709 {
00710   RooAbsArg *theVar(0);
00711   RooAbsArg *theFVar(0); // fundamental counterpart
00712   TString fullName, fullTitle, varType, myName, myTitle, Unit;
00713   rarStrParser varStrParser=getVarTNTU
00714     (varStr, option, &fullName,&fullTitle, &varType, &myName,&myTitle, &Unit);
00715   
00716   // check if it has been created
00717   theVar=getAbsVar(fullName);
00718   // now check if the short name is actually created
00719   if (!theVar) theVar=getAbsVar(myName);
00720   if (theVar) {
00721     if ("RooRealVar"==TString(theVar->ClassName())) {
00722       // add it to the params list if it is not in obs list
00723       if (_fullObs->find(theVar->GetName())) {
00724         addToObs((RooRealVar*)theVar);
00725       } else {
00726         addToParams((RooRealVar*)theVar);
00727       }
00728     }
00729     return theVar;
00730   }
00731   
00732   if ("nsRooRealVar"==varType) { // new schema
00733     theVar=new RooRealVar(myName, fullTitle, 0, Unit);
00734     // set label
00735     ((RooRealVar*)theVar)->setPlotLabel(myTitle);
00736     // construct the RRV output stream
00737     TString rrvStr=myName+"=";
00738     while (varStrParser.nArgs()>0) {
00739       rrvStr+=" "+varStrParser[0];
00740       varStrParser.Remove();
00741     }
00742     stringstream rrvStream(rrvStr.Data());
00743     RooArgSet varSet(*theVar);
00744     varSet.readFromStream((istream&)rrvStream, kFALSE);
00745     addToParams((RooRealVar*)theVar);
00746   } else if ("RooRealVar"==varType) {  // dealing with RooRealVar
00747     Double_t min(0), max(0), val(0);
00748     Int_t nBins(100);
00749     TString unit="";
00750     Bool_t setToConstant(kFALSE);
00751     if (("C"==varStrParser[0])||("c"==varStrParser[0])) {
00752       varStrParser.Remove();
00753       setToConstant=kTRUE;
00754     }
00755     if (varStrParser.nArgs()>=5) { // #1
00756       val=atof(varStrParser[0]);
00757       min=atof(varStrParser[1]);
00758       max=atof(varStrParser[2]);
00759       TString binStr=varStrParser[3](2, varStrParser[3].Length());
00760       nBins=atoi(binStr);
00761       unit=varStrParser[4];
00762     } else if (4==varStrParser.nArgs()) {
00763       if (isNumber(varStrParser[2])) {
00764         val=atof(varStrParser[0]);
00765         min=atof(varStrParser[1]);
00766         max=atof(varStrParser[2]);
00767         if (varStrParser[3].Contains("B(")) { // #2
00768           TString binStr=varStrParser[3](2, varStrParser[3].Length());
00769           nBins=atoi(binStr);
00770         } else { // #3
00771           unit=varStrParser[3];
00772         }
00773       } else { // #4
00774         min=atof(varStrParser[0]);
00775         max=atof(varStrParser[1]);
00776         TString binStr=varStrParser[2](2, varStrParser[2].Length());
00777         nBins=atoi(binStr);
00778         val=(max+min)/2.;
00779         unit=varStrParser[3];
00780       }
00781     } else if (3==varStrParser.nArgs()) {
00782       if (varStrParser[2].Contains("B(")) { // #5
00783         min=atof(varStrParser[0]);
00784         max=atof(varStrParser[1]);
00785         TString binStr=varStrParser[2](2, varStrParser[2].Length());
00786         nBins=atoi(binStr);
00787         val=(max+min)/2.;
00788       } else if (!isNumber(varStrParser[2])) { // #6
00789         min=atof(varStrParser[0]);
00790         max=atof(varStrParser[1]);
00791         val=(max+min)/2.;
00792         unit=varStrParser[2];
00793       } else { // #7
00794         val=atof(varStrParser[0]);
00795         min=atof(varStrParser[1]);
00796         max=atof(varStrParser[2]);
00797       }
00798     } else if (2==varStrParser.nArgs()) {
00799       // check if varStrParser[1] is unit
00800       if (isNumber(varStrParser[1])) { // #8
00801         min=atof(varStrParser[0]);
00802         max=atof(varStrParser[1]);
00803         val=(max+min)/2.;
00804       } else { // #9
00805         val=atof(varStrParser[0]);
00806         unit=varStrParser[1];
00807       }
00808     } else if (1==varStrParser.nArgs()) { // #10
00809       val=atof(varStrParser[0]);
00810     }
00811     if (min==max) {
00812       theVar=new RooRealVar(myName, fullTitle, val, unit);
00813     } else {
00814       theVar=new RooRealVar(myName, fullTitle, val, min, max, unit);
00815     }
00816     // set #bins
00817     ((RooRealVar*)theVar)->setBins(nBins);
00818     // set label
00819     ((RooRealVar*)theVar)->setPlotLabel(myTitle);
00820     // set error
00821     Double_t theError(0);
00822     if (theError<=0) theError=(max-min)/2.;
00823     if (theError>100) theError=100;
00824     if (min!=max) ((RooRealVar*)theVar)->setError(theError);
00825     if (setToConstant) ((RooRealVar*)theVar)->setConstant(kTRUE);
00826     // add it to the params list
00827     addToParams((RooRealVar*)theVar);
00828   } else if ("RooConstVar"==varType) {// RooConstVar
00829     Double_t val=0;
00830     if (varStrParser.nArgs()>0) val=atof(varStrParser[0]);
00831     theVar=new RooConstVar(myName, fullTitle, val);
00832   } else if ("RooUnblindOffset"==varType) {// RooUnblindOffset
00833     if (varStrParser.nArgs()<=2) {
00834       cout<<"\""<<varStr<<"\""<<endl
00835           <<"At least 3 more parameters needed"<<endl;
00836       exit(-1);
00837     }
00838     TString blindString=varStrParser[0];
00839     varStrParser.Remove();
00840     Double_t scale=atof(varStrParser[0]);
00841     varStrParser.Remove();
00842     TString blindValueName=varStrParser[0];
00843     varStrParser.Remove();
00844     RooAbsReal *blindValue=(RooAbsReal *)
00845       createAbsReal(blindValueName, blindValueName);
00846     RooAbsCategory *blindState(0);
00847     if (varStrParser.nArgs()>1) { // have blindState cat
00848       TString blindStateName=varStrParser[0];
00849       varStrParser.Remove();
00850       blindState=(RooAbsCategory *)
00851         createAbsReal(blindStateName, blindStateName);
00852     }
00853     if (blindState) {
00854       theVar=new RooUnblindOffset
00855         (myName, fullTitle, blindString, scale, *blindValue, *blindState);
00856     } else {
00857       theVar=new RooUnblindOffset
00858         (myName, fullTitle, blindString, scale, *blindValue);
00859     }
00860   } else if ("RooUnblindPrecision"==varType) {// RooUnblindPrecision
00861     if (varStrParser.nArgs()<=4) {
00862       cout<<"\""<<varStr<<"\""<<endl
00863           <<"At least 5 more parameters needed"<<endl;
00864       exit(-1);
00865     }
00866     TString blindString=varStrParser[0];
00867     varStrParser.Remove();
00868     Double_t centralValue=atof(varStrParser[0]);
00869     varStrParser.Remove();
00870     Double_t scale=atof(varStrParser[0]);
00871     varStrParser.Remove();
00872     TString blindValueName=varStrParser[0];
00873     varStrParser.Remove();
00874     RooAbsReal *blindValue=(RooAbsReal *)
00875       createAbsReal(blindValueName, blindValueName);
00876     RooAbsCategory *blindState(0);
00877     if (varStrParser.nArgs()>1) { // have blindState cat
00878       TString blindStateName=varStrParser[0];
00879       varStrParser.Remove();
00880       blindState=(RooAbsCategory *)
00881         createAbsReal(blindStateName, blindStateName);
00882     }
00883     Bool_t sin2betaMode(kFALSE);
00884     TString sin2betaModeStr=varStrParser[0];
00885     varStrParser.Remove();
00886     if (("kTRUE"==sin2betaModeStr) ||
00887         ("yes"==sin2betaModeStr) ||
00888         (0!=atoi(sin2betaModeStr))) {
00889       sin2betaMode=kTRUE;
00890     }
00891     if (blindState) {
00892       theVar=new
00893         RooUnblindPrecision(myName, fullTitle, blindString, centralValue,
00894                             scale, *blindValue, *blindState, sin2betaMode);
00895     } else {
00896       theVar=new
00897         RooUnblindPrecision(myName, fullTitle, blindString, centralValue,
00898                             scale, *blindValue, sin2betaMode);
00899     }
00900   } else if ("RooCategory"==varType) { // dealing with RooCategory
00901     if (varStrParser.nArgs()<=1) {
00902       cout<<"\""<<varStr<<"\""<<endl
00903           <<"At least 2 more parameters needed"<<endl;
00904       exit(-1);
00905     }
00906     TString idxOption=varStrParser[0];
00907     varStrParser.Remove();
00908     Bool_t useIdx(kFALSE);
00909     if (isNumber(idxOption)) {
00910       Int_t nCat=atoi(idxOption);
00911       Int_t nArgs=varStrParser.nArgs();
00912       if (nArgs==nCat*2) useIdx=kTRUE;
00913     } else {
00914       if ("noIdx"!=idxOption) useIdx=kTRUE;
00915     }
00916     if (useIdx && (varStrParser.nArgs()<=1)) {
00917       cout<<"Please give also index of cat"<<endl;
00918       exit(-1);
00919     }
00920     theVar=new RooCategory(myName, fullTitle);
00921     Int_t nTokenPerCat=1;
00922     if (useIdx) nTokenPerCat=2;
00923     Int_t nCat=varStrParser.nArgs()/nTokenPerCat;
00924     for (Int_t j=0; j<nCat; j++) {
00925       if (useIdx) { // #1
00926         ((RooCategory*)theVar)->defineType(varStrParser[j*nTokenPerCat],
00927                                          atoi(varStrParser[1+j*nTokenPerCat]));
00928       } else { // #2
00929         ((RooCategory*)theVar)->defineType(varStrParser[j*nTokenPerCat]);
00930       }
00931     }
00932     _rarCats.add(*theVar); // remember this category
00933   } else if ("RooMappedCategory"==varType) {
00934     if (varStrParser.nArgs()<=2) {
00935       cout<<"\""<<varStr<<"\""<<endl
00936           <<"At least 3 more parameters needed"<<endl;
00937       exit(-1);
00938     }
00939     TString idxOption=varStrParser[0];
00940     varStrParser.Remove();
00941     Bool_t useIdx(kFALSE);
00942     if ("noIdx"!=idxOption) useIdx=kTRUE;
00943     if (useIdx && (varStrParser.nArgs()<=2)) {
00944       cout<<"Please give default index"<<endl;
00945       exit(-1);
00946     }
00947     TString inputCatName=varStrParser[0];
00948     varStrParser.Remove();
00949     RooAbsCategory *inputCat=dynamic_cast<RooAbsCategory*>
00950       (createAbsVar(inputCatName+" "+
00951                     readConfStr(inputCatName, "", getVarSec())));
00952     if (!inputCat) {
00953       cout<<" Can not find cat named "<<inputCatName<<endl;
00954       exit(-1);
00955     }
00956     TString defCatName=varStrParser[0];
00957     varStrParser.Remove();
00958     Int_t defCatIdx=RooMappedCategory::NoCatIdx;
00959     if (useIdx) {
00960       defCatIdx=atoi(varStrParser[0]);
00961       varStrParser.Remove();
00962     }
00963     theVar=new RooMappedCategory(myName, fullTitle, *inputCat,
00964                                  defCatName, defCatIdx);
00965     Int_t nTokenPerMap=2;
00966     if (useIdx) nTokenPerMap=3;
00967     Int_t nMap=varStrParser.nArgs()/nTokenPerMap;
00968     for (Int_t j=0; j<nMap; j++) {
00969       Int_t catIdx=RooMappedCategory::NoCatIdx;
00970       if (useIdx) catIdx=atoi(varStrParser[2+j*nTokenPerMap]);
00971       ((RooMappedCategory*)theVar)->map(varStrParser[j*nTokenPerMap],
00972                                         varStrParser[1+j*nTokenPerMap],catIdx);
00973     }
00974     if (_createFundamental) {
00975       theFVar=(RooCategory*) theVar->createFundamental();
00976       _rarCats.add(*theFVar); // remember this category
00977     } else _rarCats.add(*theVar); // remember this category
00978   } else if ("RooThresholdCategory"==varType) {
00979     if (varStrParser.nArgs()<=2) {
00980       cout<<"\""<<varStr<<"\""<<endl
00981           <<"At least 3 more parameters needed"<<endl;
00982       exit(-1);
00983     }
00984     TString idxOption=varStrParser[0];
00985     varStrParser.Remove();
00986     Bool_t useIdx(kFALSE);
00987     if ("noIdx"!=idxOption) useIdx=kTRUE;
00988     if (useIdx && (varStrParser.nArgs()<=2)) {
00989       cout<<"Please give default index"<<endl;
00990       exit(-1);
00991     }
00992     TString inputVarName=varStrParser[0];
00993     varStrParser.Remove();
00994     RooAbsReal *inputVar=(RooAbsReal *)
00995       createAbsVar(inputVarName+" "+
00996                    readConfStr(inputVarName, "", getVarSec()));
00997     TString defCatName=varStrParser[0];
00998     varStrParser.Remove();
00999     Int_t defCatIdx=0;
01000     if (useIdx) {
01001       defCatIdx=atoi(varStrParser[0]);
01002       varStrParser.Remove();
01003     }
01004     theVar=new RooThresholdCategory(myName, fullTitle, *inputVar,
01005                                     defCatName, defCatIdx);
01006     Int_t nTokenPerThres=2;
01007     if (useIdx) nTokenPerThres=3;
01008     Int_t nThres=varStrParser.nArgs()/nTokenPerThres;
01009     for (Int_t j=0; j<nThres; j++) {
01010       Int_t catIdx=-99999;
01011       if (useIdx) catIdx=atoi(varStrParser[2+j*nTokenPerThres]);
01012       ((RooThresholdCategory*)theVar)->
01013         addThreshold(atof(varStrParser[j*nTokenPerThres]),
01014                      varStrParser[1+j*nTokenPerThres],catIdx);
01015     }
01016     if (_createFundamental) {
01017       theFVar=(RooCategory*) theVar->createFundamental();
01018       _rarCats.add(*theFVar); // remember this category
01019     } else _rarCats.add(*theVar); // remember this category
01020   } else if ("RooSuperCategory"==varType) {
01021     if (varStrParser.nArgs()<=1) {
01022       cout<<"\""<<varStr<<"\""<<endl
01023           <<"At least 2 more parameters needed"<<endl;
01024       exit(-1);
01025     }
01026     RooArgSet inputCats;
01027     while (varStrParser.nArgs()>0) {
01028       TString inputCatName=varStrParser[0];
01029       varStrParser.Remove();
01030       RooAbsCategory *inputCat=dynamic_cast<RooAbsCategory*>
01031         (createAbsVar(inputCatName+" "+
01032                       readConfStr(inputCatName, "", getVarSec())));
01033       if (!inputCat) {
01034         cout<<" Can not find cat named "<<inputCatName<<endl;
01035         exit(-1);
01036       }
01037       inputCats.add(*inputCat);
01038     }
01039     theVar=new RooSuperCategory(myName, fullTitle, inputCats);
01040     if (_createFundamental) {
01041       theFVar=(RooCategory*) theVar->createFundamental();
01042       _rarCats.add(*theFVar); // remember this category
01043     } else _rarCats.add(*theVar); // remember this category
01044   } else if ("RooStringVar"==varType) { // dealing with RooStringVar
01045     TString val="";
01046     if (varStrParser.nArgs()) val=varStrParser[0];
01047     theVar=new RooStringVar(myName, fullTitle, val);
01048   } else if ("RooFormulaVar"==varType) { // dealing with RooFormulaVar
01049     theVar=new RooFormulaVar(myName, myTitle, *getFormulaArgs(varStrParser));
01050     if (_createFundamental) {
01051       theFVar=(RooRealVar*) theVar->createFundamental();
01052       // check if we specify range
01053       Bool_t gotMin(kFALSE);
01054       Double_t min(0), max(0);
01055       while (varStrParser.nArgs()>0) {
01056         if (!isNumber(varStrParser[0])) {
01057           varStrParser.Remove();
01058           continue;
01059         }
01060         if (!gotMin) {
01061           min=atof(varStrParser[0]);
01062           gotMin=kTRUE;
01063         } else {
01064           max=atof(varStrParser[0]);
01065         }
01066         varStrParser.Remove();
01067       }
01068       if (gotMin) {
01069         if (min>max) {
01070           Double_t v=min;
01071           min=max;
01072           max=v;
01073         }
01074         ((RooRealVar*)theFVar)->setRange(min, max);
01075         ((RooRealVar*)theFVar)->setVal((min+max)/2.);
01076       }
01077     }
01078   } else { // not implemented for this type
01079     cout<<"\""<<varStr<<"\""<<endl
01080         <<"Var type \""<<varType<<"\" not implemented"<<endl;
01081     exit(-1);
01082   }
01083   
01084   theVar->SetName(fullName);
01085   if (theFVar) {
01086     theFVar->SetName(fullName);
01087     _rarVars.Add(theFVar);
01088     _rarOVars.Add(theVar);
01089   } else {
01090     _rarVars.Add(theVar);
01091   }
01092   
01093   cout<<varType<<" "<<fullName<<" created:\t"<<varStr<<endl;
01094   return theVar;
01095 }
01096 
01104 RooAbsArg *rarConfig::getAbsVar(TString varName)
01105 {
01106   RooAbsArg *theVar(0);
01107   if (_createFundamental) {
01108     // for _rarOVars;
01109     // first check the full name
01110     theVar=(RooAbsArg*)_rarOVars.FindObject(getFullVarName(varName));
01111     if (theVar) return theVar;
01112     // then check for varName as full name
01113     theVar=(RooAbsArg*)_rarOVars.FindObject(varName);
01114     if (theVar) return theVar;
01115   }
01116   // for _rarVars
01117   // first check the full name
01118   theVar=(RooAbsArg*)_rarVars.FindObject(getFullVarName(varName));
01119   if (theVar) return theVar;
01120   // then check for varName as full name
01121   theVar=(RooAbsArg*)_rarVars.FindObject(varName);
01122   if (theVar) return theVar;
01123   // now check if varName is actually pre-defined config
01124   // for example, mean, sigma, etc.
01125   // in this case, the AbsVar returned should be what varName refer to
01126   TString varConfigStr=readConfStr(varName, "notSet", getVarSec());
01127   if ("notSet"==varConfigStr) return theVar;
01128   TString refName=getFullVarName(varName, varName+" "+varConfigStr);
01129   if (varName==refName) return theVar; // safeguard for self refer
01130   if (refName==getFullVarName(varName)) return theVar;
01131   return getAbsVar(refName);
01132 }
01133 
01143 RooAbsArg *rarConfig::createAbsVars(TString configName,
01144                                       RooAbsCollection *argCollA,
01145                                       RooAbsCollection *argCollB)
01146 {
01147   RooAbsArg *theVar(0);
01148   // get number of vars
01149   TString varsStr=readConfStr(configName, "notSet", getVarSec());
01150   if ("notSet"==varsStr) return theVar;
01151   rarStrParser varsStrParser=varsStr;
01152   // is the first arg fullNamed
01153   Bool_t fullNamed(kFALSE);
01154   if ("fullNamed"==varsStrParser[0]) {
01155     varsStrParser.Remove();
01156     fullNamed=kTRUE;
01157   }
01158   Int_t nVar=varsStrParser.nArgs();
01159   if (nVar>0) {
01160     for (Int_t i=0; i<nVar; i++)
01161       readConfStr(varsStrParser[i], "notSet", getVarSec());
01162     cout<<"Variables defined with config \""<<configName<<"\""
01163         <<" in config file for "<<GetName();
01164     if (fullNamed) cout<<" (fullNamed)";
01165     cout<<":"<<endl;
01166     for (Int_t i=0; i<nVar; i++) {
01167       cout<<Form(" #%02d ", i)<<varsStrParser[i]<<" "
01168           <<readConfStr(varsStrParser[i], "notSet", getVarSec())<<endl;
01169     }
01170   }
01171   for (Int_t i=0; i<nVar; i++) {
01172     if (!fullNamed) theVar=createAbsReal(varsStrParser[i], varsStrParser[i]);
01173     else { // created as if the fullName is specified
01174       TString varConfigStr=
01175         readConfStr(varsStrParser[i], varsStrParser[i], getVarSec());
01176       rarStrParser varConfigStrParser=varConfigStr;
01177       TString fullName=varsStrParser[i];
01178       TString option ="N "+fullName;
01179       if (isVarType(varConfigStrParser[0])) fullName+=" "+fullName;
01180       theVar=createAbsVar(fullName+" "+varConfigStr, option);
01181     }
01182     if (argCollA) argCollA->add(*theVar);
01183     if (argCollB) argCollB->add(*theVar);
01184   }
01185   
01186   return theVar;
01187 }
01188 
01194 void rarConfig::setColLimits(RooDataSet *data, Bool_t setLimits)
01195 {
01196   if (!data) return;
01197   RooArgSet *addOns=getAddOnCols();
01198   if (!addOns) return;
01199   RooArgList addOnsList(*addOns);
01200   RooArgSet theRow(*data->get());
01201   Int_t nAddOns=addOnsList.getSize();
01202   for (Int_t i=0; i<nAddOns; i++) {
01203     TString theName=addOnsList[i].GetName();
01204     RooRealVar *theVar=dynamic_cast<RooRealVar*>(theRow.find(theName));
01205     if (!theVar) continue;
01206     RooRealVar *fVar=dynamic_cast<RooRealVar*>(_rarVars.FindObject(theName));
01207     if (!fVar) {
01208       cout<<" No fundamental var created for "<<theName<<endl;
01209       exit(-1);
01210     }
01211     if (setLimits) {
01212       Double_t theMin=fVar->getMin();
01213       Double_t theMax=fVar->getMax();
01214       Double_t dMin, dMax;
01215       data->getRange(*theVar, dMin, dMax);
01216       if (dMin<theMin) {
01217         cout<<" W A R N I N G ! ! !"<<endl
01218             <<" The allowable min value ("<<theMin<<") for "<<theVar->GetName()
01219             <<" is larger than the min value ("<<dMin<<")"
01220             <<" in dataset "<<data->GetName()<<endl
01221             <<" The min of "<<theVar->GetName()<<" is set to "<<dMin
01222             <<" for "<<data->GetName()<<endl
01223             <<" Please make sure the range for "<<theVar->GetName()
01224             <<" is valid"<<endl
01225             <<" Your results based on "<<data->GetName()
01226             <<" may NOT be accurate!"<<endl
01227             <<endl;
01228         theMin=dMin;
01229       }
01230       if (dMax>theMax) {
01231         cout<<" W A R N I N G ! ! !"<<endl
01232             <<" The allowable max value ("<<theMax<<") for "<<theVar->GetName()
01233             <<" is smaller than the max value ("<<dMax<<")"
01234             <<" in dataset "<<data->GetName()<<endl
01235             <<" The max of "<<theVar->GetName()<<" is set to "<<dMax
01236             <<" for "<<data->GetName()<<endl
01237             <<" Please make sure the range for "<<theVar->GetName()
01238             <<" is valid"<<endl
01239             <<" Your results based on "<<data->GetName()
01240             <<" may NOT be accurate!"<<endl
01241             <<endl;
01242         theMax=dMax;
01243       }
01244       theVar->setRange(theMin, theMax);
01245     } else
01246       theVar->removeRange();
01247     //cout<<" Derived obs "<<theName <<" set to ("<<theVar->getMin()<<", "
01248     //<<theVar->getMax()<<")"<<endl;
01249   }
01250 }
01251 
01259 void rarConfig::addColumns(RooDataSet *data, Bool_t addColmns,
01260                            Bool_t setLimits)
01261 {
01262   if (!data) return;
01263   RooArgSet *addOns=getAddOnCols();
01264   if (!addOns) return;
01265   // add the columns
01266   if (addColmns) {
01267     // data->addColumns(*addOns); // gives errors "can only add to owned list"
01268     TIterator* iter = addOns->createIterator();
01269     RooAbsArg *temp(0);
01270     while ( temp = (RooAbsArg*)(iter->Next()) ) {
01271       data->addColumn(*temp);
01272     }
01273   }
01274   // then set the limits for those obs if it is RRV
01275   setColLimits(data, setLimits);
01276 }
01277 
01293 RooDataSet *rarConfig::createDataSet(const char *dsStr, Bool_t &isUB, TString wgtVarName)
01294 {
01295   isUB=kFALSE;
01296   rarStrParser datasetStrParser=dsStr;
01297   TString myName=datasetStrParser[0];
01298   datasetStrParser.Remove();
01299   TString dsType=datasetStrParser[0];
01300   datasetStrParser.Remove();
01301   TString myTitle=datasetStrParser[0];
01302   datasetStrParser.Remove();
01303   
01304   // cout<<endl<<" For "<<myName<<endl;
01305   RooDataSet *data(0);
01306   if ("ascii"==dsType) { // dealing with ascii text file method
01307     RooDataSet *theData(0); // temporary data set
01308     if (1==datasetStrParser.nArgs()) { // #1
01309       theData=RooDataSet::read(datasetStrParser[0], *getPrimaryObs());
01310     } else if(2==datasetStrParser.nArgs()) { // #2, with options
01311       theData=RooDataSet::read(datasetStrParser[0], *getPrimaryObs(),
01312                             datasetStrParser[1]);
01313     } else if(3==datasetStrParser.nArgs()) { // #3, with multi text files
01314       theData=RooDataSet::read(datasetStrParser[0], *getPrimaryObs(),
01315                             datasetStrParser[1], datasetStrParser[2]);
01316     } else if(4==datasetStrParser.nArgs()) { // #4, with indexCatName
01317       theData=RooDataSet::read(datasetStrParser[0], *getPrimaryObs(),
01318                             datasetStrParser[1], datasetStrParser[2],
01319                             datasetStrParser[3]);
01320     } else {
01321       cout<<"Not implemented for "<<datasetStrParser.nArgs()<<" args"<<endl;
01322       exit(-1);
01323     }
01324     // if weight
01325     // RooArgSet primaryObs = *getPrimaryObs();
01326     // primaryObs.add(evtWgt); // add the weight to the list
01327     // data = new RooDataSet("theData", "theData", tree, primaryObs, cut, wgtVarName);
01328     // TFile f("/tmp/temporay.root");
01329     // data->Write();
01330     // f.ls();
01331     // f.Close();
01332     // RooDataSet *d = (RooDataSet*) f.FindObject("d");
01333     addColumns(theData);
01334     // copy and add weight variable
01335     data=new RooDataSet("asciiData", "asciiData", theData, *_fullObs, "1", wgtVarName);
01336   } else if ("root"==dsType) { // dealing with root file
01337     if (datasetStrParser.nArgs()<2) {
01338       cout<<"At least two more args needed"<<endl;
01339       exit(-1);
01340     }
01341     TString cut="";
01342     if (datasetStrParser.nArgs()>2) cut=datasetStrParser[2];
01343  
01344     // open the file
01345     TFile f(datasetStrParser[0]);
01346 
01347     TTree *tree = (TTree*) f.Get(datasetStrParser[1]);
01348     if (!tree) {
01349       cout<< "Can not find tree "<<datasetStrParser[1] <<" in file "<<datasetStrParser[0]<<endl;
01350       exit(-1);
01351     }
01352     
01353     // now create the dataset and add weight variable at same time
01354     RooArgSet primaryObs = *getPrimaryObs();
01355     if (wgtVarName != "") {
01356       RooRealVar evtWgt(wgtVarName, "weight", 1.0);
01357       primaryObs.add(evtWgt); // add the weight to the list
01358       data = new RooDataSet("theData", "theData", tree, primaryObs, cut, wgtVarName);
01359     } else {
01360       data = new RooDataSet("theData", "theData", tree, primaryObs, cut);
01361     }
01362 
01363     addColumns(data);
01364     //data->Print("v");
01365   } else if ("add"==dsType) { // dealing with add method
01366     // create dataset by adding existing datasets
01367     cout<<"RooDataSet::add: adding";
01368     for (Int_t i=0; i<datasetStrParser.nArgs(); i++) {
01369       cout<<" "<<datasetStrParser[i];
01370     }
01371     cout<<endl;
01372     data=new RooDataSet("addData", "add dataset", *_fullObs);
01373     Int_t nDS=datasetStrParser.nArgs()/2;
01374     for (Int_t i=0; i<nDS; i++) {
01375       // get dataset #i
01376       RooDataSet *theData=getData(datasetStrParser[0]);
01377       datasetStrParser.Remove();
01378       if (!theData) continue;
01379       Int_t srcNevt=(Int_t)theData->numEntries();
01380       Double_t nEvtScale=atof(datasetStrParser[0]);
01381       datasetStrParser.Remove();
01382       Int_t nEvt=(Int_t)nEvtScale;
01383       if (nEvtScale<1) { // if less than 1, it is scaling factor
01384         nEvt=(Int_t)(.5+nEvtScale*srcNevt);
01385       }
01386       if ((nEvt<=0) || (nEvt>srcNevt)) nEvt=srcNevt;
01387       if ((nEvt>0)&&(nEvt<srcNevt)) isUB=kTRUE;
01388       cout<<" Adding "<<nEvt<<" events from "<<theData->GetName()<<endl;
01389       // create index vector
01390       vector<Int_t> indexVector;
01391       for (Int_t j=0; j<srcNevt; j++) {
01392         indexVector.push_back(j);
01393       }
01394       for (Int_t j=0; j<nEvt; j++) {
01395         // random access to index to make sure the sample added
01396         // is not sequence dependent to the src data
01397         Int_t randomIdx=
01398           RooRandom::randomGenerator()->Integer(indexVector.size());
01399         RooArgSet *theRow=(RooArgSet *)theData->get(indexVector[randomIdx]);
01400         data->add(*theRow);
01401         // remove the index so it won't be selected again
01402         vector<Int_t>::iterator the_iterator;
01403         the_iterator=indexVector.begin();
01404         the_iterator+=randomIdx;
01405         indexVector.erase(the_iterator);
01406       }
01407     }
01408     addColumns(data);
01409     cout<<"RooDataSet::add: added "<<data->numEntries()<<" events"<<endl;
01410   } else if ("reduce"==dsType) { // dealing with reduce method
01411     if (datasetStrParser.nArgs()<=0) {
01412       cout<<"need more args"<<endl;
01413       exit(-1);
01414     }
01415     TString srcDataName=datasetStrParser[0];
01416     datasetStrParser.Remove();
01417     RooDataSet *theData=getData(srcDataName);
01418     if (!theData) {
01419       cout<<"Can not find dataset named \""<<srcDataName<<"\""<<endl;
01420       exit(-1);
01421     }
01422     TString cuts="1";
01423     if (datasetStrParser.nArgs()>=1) {
01424       cuts=datasetStrParser[0];
01425       datasetStrParser.Remove();
01426     }
01427     data=(RooDataSet*)theData->reduce(cuts);
01428 
01429     cout<<"RooDataSet::reduce: reduced "<<data->numEntries()<<" events"<<endl;
01430  
01431   } else if ("hist"==dsType) { // dealing with RooDataHist
01432     if (datasetStrParser.nArgs()<=0) {
01433       cout<<"need more args"<<endl;
01434       exit(-1);
01435     }
01436     TString srcDataName=datasetStrParser[0];
01437     datasetStrParser.Remove();
01438     RooDataSet *theData=getData(srcDataName);
01439     if (!theData) {
01440       cout<<"Can not find dataset named \""<<srcDataName<<"\""<<endl;
01441       exit(-1);
01442     }
01443     data=(RooDataSet*)
01444       new RooDataHist("theData", "theData", *getPrimaryObs(), *theData);
01445   } else {
01446     cout<<"Not implemented for "<<datasetStrParser[1]<<" dataset type"<<endl;
01447     exit(-1);
01448   }
01449   
01450   // change its name and title
01451   data->SetName(myName);
01452   data->SetTitle(myTitle);
01453   
01454   return data;
01455 }
01456 
01462 void rarConfig::computeCorrelations(RooArgList varList,
01463                                     const RooDataSet *data) {
01464   RooArgSet theFVars(*data->get());
01465   // remove any categories
01466   RooArgSet theVars;
01467   for (Int_t i=0; i<varList.getSize(); i++) {
01468     RooAbsArg *theVar=varList.at(i);
01469     if (_rarOVars.FindObject(theVar->GetName())) continue;
01470     if (getCats()->find(theVar->GetName())) continue;
01471     theVars.add(*(theFVars.find(theVar->GetName())));
01472   }
01473   // make sure there are vars for correlation matrix
01474   if (theVars.getSize()<=1) return;
01475   
01476   // the following codes are mainly from Q2BTimeMLFit::computeCorrelations
01477   unsigned i,j,k;
01478   unsigned nEvt = data->numEntries();
01479   unsigned nPar =  theVars.getSize();
01480   double* average = new double[nPar];
01481   RooAbsReal **arg = new RooAbsReal*[nPar];
01482   TIterator* iter = theVars.createIterator();
01483   i=0;
01484   RooAbsArg *y(0);
01485   while ( (y=(RooAbsArg*)(iter->Next())) ) {
01486     arg[i] =  dynamic_cast<RooAbsReal*>(y);
01487     average[i]=0;
01488     i++;
01489   }
01490   for (i=0;i<nEvt;++i) {
01491     data->get(i);
01492     for (unsigned j=0;j<nPar;++j)
01493       if (arg[j]!=0) average[j]+=(arg[j]->getVal())/nEvt;
01494   }
01495   double** cor = new double*[nPar];
01496   for (i=0;i<nPar;++i) {
01497     cor[i] = new double[nPar];
01498     for (unsigned j=0;j<nPar;++j) cor[i][j]=0;
01499   }
01500   for (i=0;i<nEvt;++i) {
01501     data->get(i);
01502     for (j=0;j<nPar;++j) for (k=0;k<nPar;++k) {
01503       if (arg[j]!=0)
01504         cor[j][k]+=(arg[j]->getVal()-average[j])*(arg[k]->getVal()-average[k]);
01505     }
01506   }
01507   for (j=0;j<nPar;++j) for (k=0;k<nPar;++k)  {
01508     if (j!=k && arg[j]!=0 && arg[k]!=0)
01509       cor[j][k]/=sqrt(cor[j][j]*cor[k][k]);
01510   }
01511   
01512   cout << endl << " Correlation matrix for " << data->GetName()
01513        << ":" << endl;
01514   cout << "        ";
01515   char chrbuf[9], label[9];
01516   for (i=0; i< nPar-1;i++) {
01517     strncpy(label, arg[i]->GetName(), 8);  label[8] = '\0';
01518     sprintf(chrbuf, "%8s", label);
01519     cout << "  " << chrbuf;
01520   }
01521   cout << endl;
01522   for (i=0;i<nPar;++i) {
01523     if (i>0) {
01524       strncpy(label, arg[i]->GetName(), 8);  label[8] = '\0';
01525       sprintf(chrbuf, "%8s", label);
01526       cout << chrbuf;
01527       for (j=0;j<i;++j) {
01528         sprintf(chrbuf, "%8.4f", cor[i][j]);
01529         cout << "  " << chrbuf;
01530       }
01531       cout << endl;
01532     }
01533     delete[] cor[i];
01534   }
01535   cout << endl;
01536   delete[] cor;
01537   delete[] average;
01538   delete[] arg;
01539 
01540   // FFW but only outside babar
01541   /*
01542   const TMatrixDSym* corn = data->correlationMatrix(theVars) ;
01543   //const TMatrixDSym* covn = data->covarianceMatrix(theVars) ;
01544 
01545   // Print correlation, covariance matrix
01546   cout << "correlation matrix for " << data->GetName() << endl ;
01547   corn.Print() ;
01548   //cout << "covariance matrix" << endl ;
01549   //covn.Print() ;
01550   */
01551 }
01552 
01576 rarBasePdf *rarConfig::createPdf(const char *configStr)
01577 {
01578   rarStrParser configStrParser=configStr;
01579   if(configStrParser.nArgs()<1) {
01580     cout<<"\""<<configStr<<"\""<<endl
01581         <<" At least 1 arg needed for rarConfig::createPdf"<<endl;
01582     exit(-1);
01583   }
01584   rarBasePdf *thePdf(0);
01585   // get the name of the pdf
01586   TString theName=configStrParser[0];
01587   configStrParser.Remove();
01588   // check if the pdf has already been created
01589   if (thePdf=(rarBasePdf*)_rarPdfs.FindObject(theName)) return thePdf;
01590   TString theConfigSec=theName+" Config";
01591   // get the pdf type
01592   if(configStrParser.nArgs()<1) {
01593     cout<<" You have to specify the pdf type"<<endl;
01594     exit(-1);
01595   }
01596   TString thePdfType=configStrParser[0];
01597   configStrParser.Remove();
01598   // get the pdf title
01599   TString theTitle(thePdfType+" "+GetTitle()); // the default title
01600   if(configStrParser.nArgs()>0) {
01601     theTitle=configStrParser[0]; // this is the title
01602     configStrParser.Remove();
01603   }
01604   const char *dummy("");
01605   RooDataSet *_theData=getData(dummy);
01606   rarDatasets *_datasets=getDatasets();
01607   // create different types of PDFs
01608   if ("MLPdf"==thePdfType) {
01609     if (getPdfType()!="MLFitter") {
01610       cout<<"Only the final mlFitter can create "<<thePdfType<<endl;
01611       exit(-1);
01612     }
01613     thePdf=new rarMLPdf(_configFile, theConfigSec, configStr,
01614                         _datasets, _theData, theName, theTitle);
01615   } else if ("ProdPdf"==thePdfType) {
01616     thePdf=new rarProd(_configFile, theConfigSec, configStr,
01617                        _datasets, _theData, theName, theTitle);
01618   } else if (("AddPdf"==thePdfType)||("AddModel"==thePdfType)) {
01619     thePdf=new rarAdd(_configFile, theConfigSec, configStr,
01620                       _datasets, _theData, theName, theTitle);
01621   } else if ("Simultaneous"==thePdfType) {
01622     thePdf=new rarSimPdf(_configFile, theConfigSec, configStr,
01623                          _datasets, _theData, theName, theTitle);
01624   } else if ("Exp"==thePdfType) {
01625     thePdf=new rarExp(_configFile, theConfigSec, configStr,
01626                       _datasets, _theData, theName, theTitle);
01627   } else if (("Gaussian"==thePdfType)||("BreitWigner"==thePdfType) ||
01628              ("Landau"==thePdfType)) {
01629     thePdf=new rarGaussian(_configFile, theConfigSec, configStr,
01630                            _datasets, _theData, theName, theTitle);
01631   } else if ("TwoGauss"==thePdfType) {
01632     thePdf=new rarTwoGauss(_configFile, theConfigSec, configStr,
01633                            _datasets, _theData, theName, theTitle);
01634   } else if ("OsipDisc"==thePdfType) {
01635     thePdf=new rarOsipDisc(_configFile, theConfigSec, configStr,
01636                            _datasets, _theData, theName, theTitle);
01637   } else if (("TriGauss"==thePdfType)||
01638              ("TriGaussModel"==thePdfType)||("GexpShape"==thePdfType)) {
01639     thePdf=new rarTriGauss(_configFile, theConfigSec, configStr,
01640                            _datasets, _theData, theName, theTitle);
01641   } else if (("BifurGauss"==thePdfType)||("BGGauss"==thePdfType)) {
01642     thePdf=new rarBifurGauss(_configFile, theConfigSec, configStr,
01643                              _datasets, _theData, theName, theTitle);
01644   } else if ("Novosibirsk"==thePdfType) {
01645     thePdf=new rarNovosibirsk(_configFile, theConfigSec, configStr,
01646                               _datasets, _theData, theName, theTitle);
01647   } else if ("Ballack"==thePdfType) {
01648     thePdf=new rarBallack(_configFile, theConfigSec, configStr,
01649                           _datasets, _theData, theName, theTitle);
01650   } else if ("CBShape"==thePdfType) {
01651     thePdf=new rarCBShape(_configFile, theConfigSec, configStr,
01652                           _datasets, _theData, theName, theTitle);
01653   } else if (("Polynomial"==thePdfType)||("Chebychev"==thePdfType)) {
01654     thePdf=new rarPoly(_configFile, theConfigSec, configStr,
01655                        _datasets, _theData, theName, theTitle);
01656   } else if ("ArgusBG"==thePdfType) {
01657     thePdf=new rarArgusBG(_configFile, theConfigSec, configStr,
01658                           _datasets, _theData, theName, theTitle);
01659   } else if ("Step"==thePdfType) {
01660     thePdf=new rarStep(_configFile, theConfigSec, configStr,
01661                        _datasets, _theData, theName, theTitle);
01662   } else if ("Binned"==thePdfType) {
01663     thePdf=new rarBinned(_configFile, theConfigSec, configStr,
01664                          _datasets, _theData, theName, theTitle);
01665   } else if (("Keys"==thePdfType)||("2DKeys"==thePdfType)) {
01666     thePdf=new rarKeys(_configFile, theConfigSec, configStr,
01667                        _datasets, _theData, theName, theTitle);
01668   } else if ("Generic"==thePdfType) {
01669     thePdf=new rarGeneric(_configFile, theConfigSec, configStr,
01670                           _datasets, _theData, theName, theTitle);
01671   } else if ("HistPdf"==thePdfType) {
01672     thePdf=new rarHistPdf(_configFile, theConfigSec, configStr,
01673                           _datasets, _theData, theName, theTitle);
01674   } else if ("GaussModel"==thePdfType) {
01675     thePdf=new rarGaussModel(_configFile, theConfigSec, configStr,
01676                              _datasets, _theData, theName, theTitle);
01677   } else if (("BCPGenDecay"==thePdfType)||("BDecay"==thePdfType)||
01678              ("Decay"==thePdfType)) {
01679     thePdf=new rarDecay(_configFile, theConfigSec, configStr,
01680                         _datasets, _theData, theName, theTitle);
01681   } else if ("Cruijff"==thePdfType) {
01682     thePdf=new rarCruijff(_configFile, theConfigSec, configStr,
01683                           _datasets, _theData, theName, theTitle);
01684   } else if ("MultPdf"==thePdfType) {
01685     thePdf=new rarMultPdf(_configFile, theConfigSec, configStr,
01686                          _datasets, _theData, theName, theTitle);
01687   } else if ("UsrPdf"==thePdfType) {
01688     thePdf=new rarUsrPdf(_configFile, theConfigSec, configStr,
01689                          _datasets, _theData, theName, theTitle);
01690   } else if ("RelBreitWigner"==thePdfType) {
01691     thePdf=new rarRelBreitWigner(_configFile, theConfigSec, configStr,
01692                          _datasets, _theData, theName, theTitle);
01693   } else if ("Lass"==thePdfType) {
01694     thePdf=new rarLass(_configFile, theConfigSec, configStr,
01695                        _datasets, _theData, theName, theTitle);
01696   } else if ("Voigtian"==thePdfType) {
01697     thePdf=new rarVoigtian(_configFile, theConfigSec, configStr,
01698                        _datasets, _theData, theName, theTitle);
01699   } else if ("GounarisSakurai"==thePdfType) {
01700     thePdf=new rarGounarisSakurai(_configFile, theConfigSec, configStr,
01701                                   _datasets, _theData, theName, theTitle);
01702   } else if ("Flatte"==thePdfType) {
01703     thePdf=new rarFlatte(_configFile, theConfigSec, configStr,
01704                          _datasets, _theData, theName, theTitle);
01705   } else if ("Uniform"==thePdfType) {
01706     thePdf=new rarUniform(_configFile, theConfigSec, configStr,
01707                          _datasets, _theData, theName, theTitle);
01708   } else if ("Threshold"==thePdfType) {
01709     thePdf=new rarThreshold(_configFile, theConfigSec, configStr,
01710                          _datasets, _theData, theName, theTitle);
01711   } else {
01712     cout<<"\""<<configStr<<"\""<<endl
01713         <<"Pdf \""<<thePdfType<<"\" not implemented"<<endl;
01714     exit(-1);
01715   }
01716   _rarPdfs.Add(thePdf);
01717   
01718   cout<<thePdfType<<" created"<<endl<<endl;
01719   return thePdf;
01720 }