// File: $Id: cpmix.cc,v 1.4 2002/02/08 23:13:59 verkerke Exp $
//
// Driver macro for sin2beta analysis
//
// Author: Wouter Verkerke
//
// (C) 2002 University of California
//


// Global variables to allow access to datasets, 
// pdf's and fitresults after running
RooDataSet *brecoRawData, *charmRawData, *allData,
           *goldData, *bmixData, *cpfitData ;
RooAbsPdf *pdfMbsB, *pdfMbsG, *pdfMIX ;
RooCPMixFit *cpfit ;
RooFitResult* r ;


void cpmix(const char* configFile="cpmix.config", const char* driverOpt="F") 
{
  // Driver options
  //   "I" - Instantiate RooCPMixFit
  //   "D" - Load datasets
  //   "P" - Build PDFs
  //   "F" - Fit PDFs to data

  cout << "cpmix: configFile = '" << configFile << "' driverOpt = " << driverOpt << endl ;
  TString dopt(driverOpt) ;
  dopt.ToUpper() ;
  RooRealVar::printScientific(1) ;

  // Setup blinding configuration
  RooStringVar blindingString("blindingString","Blinding String","") ;
  RooCategory blindingStatus("blindingStatus","Blinding Status") ;
  RooRealVar blindingSigma("blindingSigma","Blinding Sigma",0) ;
  RooRealVar blindingMPoint("blindingMPoint","Blinding Mirror Point",0) ;
  blindingStatus.defineType("Blind") ;
  blindingStatus.defineType("Unblind") ;
  RooArgSet blindConfig(blindingString,blindingStatus,blindingSigma,blindingMPoint) ;
  blindConfig.readFromFile(configFile,"READ","Blinding") ;
  RooArgSet* blindNotRead = blindConfig.selectByAttrib("READ",kFALSE) ;
  if (blindNotRead->getSize()>0) {
    cout << "ERROR: Blinding configuration not completely read, missing parts are " ;
    blindNotRead->Print("1") ;
    return ;
  }

  // Instantiate RooCPMixFit object (owns all sin2beta PDF/dataset prototypes)
  cpfit = new RooCPMixFit(blindingString.getVal(),blindingMPoint.getVal(),blindingSigma.getVal()) ;

  if (dopt.Contains("I")) return ;

  // Instantiate PDF builder with sin2beta prototypes
  RooSimPdfBuilder* mgr = new RooSimPdfBuilder(cpfit.protoPdfSet()) ;


  // Read the [Input Data] section of the configuration file
  RooCategory inputSource("inputSource","input source") ;
  inputSource.defineType("ROOT") ;
  inputSource.defineType("ASCII") ;
  RooStringVar rootFile("rootFile","ROOT file with persisted allData dataset","") ;
  RooStringVar charmDir("charmDir","Directory with Charmonium ASCII files","") ;
  RooStringVar charmFiles("charmFiles","List of Charmonium ASCII files","",2048) ;
  RooStringVar brecoDir("brecoDir","Directory with BReco ASCII files","") ;
  RooStringVar brecoFiles("brecoFiles","List of BReco ASCII files","",2048) ;
  RooStringVar dsel_prefitGold("dsel_prefitGold","Data selection for CP golden modes prefit","1") ;
  RooStringVar dsel_prefitBMix("dsel_prefitBMix","Data selection for Mixing modes prefit","1") ;
  RooStringVar dsel_cpfit("dsel_cpfit","Data selection for CP/Mixing fit","1") ;

  RooArgSet fileConfig(inputSource,rootFile,charmDir,charmFiles,brecoDir,brecoFiles,dsel_prefitGold,dsel_prefitBMix,dsel_cpfit) ;
  fileConfig.readFromFile(configFile,0,"Input Data") ;

  // Root or ASCII input?
  if (inputSource=="ASCII") {
    // Read the ASCII files listed in [Input Data]
    cout << "cpmix: reading events from ASCII source" << endl ;

    // Make the raw datasets from the ascii files
    charmRawData = cpfit.loadAsciiFiles(charmFiles.getVal(),charmDir.getVal(),"q") ;
    brecoRawData = cpfit.loadAsciiFiles(brecoFiles.getVal(),brecoDir.getVal(),"q") ;  
    
    // Make a combined dataset of selected ascii variables
    allData = (RooDataSet*) charmRawData->reduce(cpfit.selDataVars) ;
    allData->append(*brecoRawData) ;
  } else {
    // Read the ROOT file listed in [Input Data]
    cout << "cpmix: reading events from ROOT source" << endl ;
    
    TFile *f = new TFile(rootFile.getVal()) ;
    allData = f->Get("allData") ;
  }

  // Make data subsets for the fits
  goldData  = (RooDataSet*) allData->reduce(dsel_prefitGold.getVal()) ;
  bmixData  = (RooDataSet*) allData->reduce(dsel_prefitBMix.getVal()) ;
  cpfitData = (RooDataSet*) allData->reduce(dsel_cpfit.getVal()) ;
  
  if (dopt.Contains("D")) return ;

  // Set delta-t range to +/- 50 emulating +/- inf  for CPExtract/tFit compatibility
  cpfit->dtReco->setFitRange(-50,50) ;
  ((RooRealVar*)goldData->get()->find("dtReco"))->setFitRange(-50,50) ;
  ((RooRealVar*)bmixData->get()->find("dtReco"))->setFitRange(-50,50) ;
  ((RooRealVar*)cpfitData->get()->find("dtReco"))->setFitRange(-50,50) ;

  // Read the PDF configurations for the mB(breco) mB(gold) and cp fit
  // from the configuration file and build the PDFs accordingly
  RooArgSet* rsbConfig = mgr->createProtoBuildConfig() ;
  rsbConfig->readFromFile(configFile,0,"mB-prefit(BMix) structure") ;
  pdfMbsB = mgr->buildPdf(*rsbConfig,cpfitData) ;			         
  rsbConfig->readFromFile(configFile,0,"mB-prefit(Gold) structure") ;
  pdfMbsG = mgr->buildPdf(*rsbConfig,cpfitData) ;
  rsbConfig->readFromFile(configFile,0,"CP/mixing-fit structure") ;
  pdfMIX  = mgr->buildPdf(*rsbConfig,cpfitData,*cpfit.tagCatL) ;

  // Read the initial parameters from 
  // configuration file section [Initial parameter values]
  RooArgSet* mixParams = pdfMIX->getParameters(cpfitData) ;
  mixParams->readFromFile(configFile,"READ","Initial parameter values") ;

  // List all paramaters as read from the file
  cout << "Initial Parameters from '" << configFile << "' : " ; 
  mixParams->selectByAttrib("READ",kTRUE)->Print("v") ;

  // List all parameters that have not been configured from file
  cout << "Parameters NOT configured from '" << configFile << "' : " ;
  mixParams->selectByAttrib("READ",kFALSE)->Print("v") ;

  if (dopt.Contains("P")) return ;

  // Run the mBshape pre-fits
  pdfMbsB->fitTo(*bmixData,"mht") ;
  pdfMbsG->fitTo(*goldData,"mht") ;

  // Fix the mBshape parameters
  mbsbVars = pdfMbsB->getParameters(allData) ;
  mbsbVars->Print("v") ;
  mbsbVars->setAttribAll("Constant",kTRUE) ;

  mbsgVars = pdfMbsG->getParameters(allData) ;
  mbsgVars->Print("v") ;
  mbsgVars->setAttribAll("Constant",kTRUE) ;

  // Optionally set all parameters except sin2beta constant
  RooCategory constSwitch("setAllParamsConstantExceptSin2b","const switch") ;
  constSwitch.defineType("Yes") ;
  constSwitch.defineType("No") ;
  constSwitch = "No" ;
  RooArgSet(constSwitch).readFromFile(configFile,"READ","Driver options") ;
  if (constSwitch=="Yes") {
    cout << "cpmix: Setting all parameters except sin2beta constant" << endl ;
    RooArgSet* mixVars = pdfMIX->getParameters(cpfitData) ;
    mixVars.remove(*mixVars->find("sin2b")) ;
    mixVars->setAttribAll("Constant",kTRUE) ;
  }

  // Perform the mixing fit
  cout << "starting cp/mixing fit" << endl ;
  r = pdfMIX->fitTo(*cpfitData,*cpfit.dtErr,"mhrvt0","cpsd") ;

  // Print out the fit result and the sin2beta correlation matrix row
  r->Print("v") ;
  (r->correlation(*cpfit.sin2b))->Print("v") ;

  return ;
}