GOmegaPiApp.cc

/**
 * @file GOmegaPiApp.cc
 */

/************************************************************************************************/
// Boost header files go here
#include "boost/date_time/posix_time/posix_time.hpp"
#include "boost/logic/tribool.hpp"

// GenEvA header files go here
#include <courtier/GAsioHelperFunctions.hpp>
#include <courtier/GAsioTCPClientT.hpp>
#include <courtier/GAsioTCPConsumerT.hpp>
#include <geneva/GBrokerEA.hpp>
#include <geneva/GEvolutionaryAlgorithm.hpp>
#include <geneva/GIndividual.hpp>
#include <geneva/GMultiThreadedEA.hpp>

#include "PwaUtils/pbarpStates.hh"
#include "Examples/pbarpToOmegaPi/AbsOmegaPiLh.hh"
#include "Examples/pbarpToOmegaPi/OmegaPiLhGamma.hh"
#include "Examples/pbarpToOmegaPi/OmegaPiLhGammaBw.hh"
//#include "Examples/pbarpToOmegaPi/OmegaPiLhOmega.hh"

#include "Examples/pbarpToOmegaPi/pbarpToOmegaPi0States.hh"
// The individual that should be optimized
#include "Examples/pbarpToOmegaPi/GOmegaPiIndividual.hh"

// Declares a function to parse the command line
#include "Examples/pbarpToOmegaPi/OmegaPiEventList.hh"
#include "Examples/pbarpToOmegaPi/OmegaPiHist.hh"

#include "Examples/pbarpToOmegaPi/OmegaPiData.hh"

#include "Setup/PwaEnv.hh"
#include "Particle/ParticleTable.hh"
#include "Particle/Particle.hh"
#include "Event/EventList.hh"
#include "Event/Event.hh"
#include "Event/CBElsaReader.hh"
#include "Particle/PdtParser.hh"

#include "ErrLogger/ErrLogger.hh"
#include "Examples/pbarpToOmegaPi/GArgumentParser.hh"

#include "Examples/pbarpToOmegaPi/MOmegaPiFcn.hh"

#include "qft++/topincludes/tensor.hh"

#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnStrategy.h"

#include <complex>

#include "Examples/pbarpToOmegaPi/spindensityhist.hh"

/************************************************************************************************/
// Namespace aliases, so we do not need to quote the entire namespace name
namespace gg = Gem::Geneva;
namespace gp = Gem::Pawian;
namespace gc = Gem::Courtier;
namespace rm = ROOT::Minuit2;
namespace bp = boost::posix_time;
namespace bl = boost::logic;

/************************************************************************************************/
/**
 * This function constructs the path to the file.
 */
void constructPath(const std::string &thePrefix, const unsigned pbarMom, std::string &outFilePath)
{
  std::stringstream sstrDatFile; //String Stream for die construction of the path to the parameter File;
  sstrDatFile << thePrefix; 
  sstrDatFile.width(4);
  sstrDatFile.fill('0');
  sstrDatFile << right << pbarMom << ".dat";
  outFilePath = sstrDatFile.str();
}


/************************************************************************************************/
/**
 * This function checks if the file in the path theFilePath exists
 */
bool checkFileExist(const std::string &theFilePath)
{
  ifstream datChk(theFilePath.c_str());
  if (datChk) { return true; } 
  else { return false; }
}

/************************************************************************************************/
/**
 * The actual optimization procedure using the Geneva library collection. Note that this function
 * will return a boost::logic::indeterminate value if this is a Geneva client in networked mode
 *
 * @return A boost::logic::tribool indicating whether the fit was successful or this is a client in networked mode
 */
bl::tribool GenEvA(
      boost::shared_ptr<const OmegaPiEventList> &theOmegaPiEventPtr
      , boost::shared_ptr<const pbarpToOmegaPi0States> &pbarpToOmegaPi0StatesPtr
      , ApplicationParameter &theAppParams
      , OmegaPiData::fitParamVal &finalFitParm
      , boost::shared_ptr<AbsOmegaPiLh> &finalOmegaPiLh
) {
  //--------------------------------------------------------------------------------------------
  // First check whether this is a client in networked mode. We then just start the listener and
  // return when it has finished.
  if(theAppParams.getParallelizationMode()==2 && !theAppParams.getServerMode()) {
    // Create a model that holds the static data needed for processing
    boost::shared_ptr<gp::GOmegaPiIndividual> gopi_ptr( new gp::GOmegaPiIndividual(finalOmegaPiLh) );

//     boost::shared_ptr<gc::GAsioTCPClientT<gg::GIndividual> > p(
// 	  new gc::GAsioTCPClientT<gg::GIndividual>(
//              theAppParams.getIp()
//              , boost::lexical_cast<std::string>(theAppParams.getPort())
// 	     , gopi_ptr
//           )
//     );
//     p->setMaxStalls(0); // An infinite number of stalled data retrievals
//     p->setMaxConnectionAttempts(100); // Up to 100 failed connection attempts

//     // Prevent return of unsuccessful mutation attempts to the server
//     p->returnResultIfUnsuccessful(theAppParams.getReturnRegardless());

//     // Start the actual processing loop
//     p->run();

    return bl::indeterminate;
  }

  //--------------------------------------------------------------------------------------------
  Info << "GenEvA fit start.\n" << endmsg;

  // Create the first set of parent individuals. Initialization of parameters is done randomly.
  std::vector<boost::shared_ptr<gp::GOmegaPiIndividual> > parentIndividuals;
  for(std::size_t p = 0 ; p<theAppParams.getNParents(); p++) {
    boost::shared_ptr<gp::GOmegaPiIndividual> gopi_ptr( new gp::GOmegaPiIndividual(finalOmegaPiLh) );
    gopi_ptr->setProcessingCycles(theAppParams.getProcessingCycles());
    // Make sure we start with a unique individual
    gopi_ptr->randomInit();
    
    parentIndividuals.push_back(gopi_ptr);
  }

  //----------------------------------------------------------------------------------------------
  // We can now start creating populations. We refer to them through the base class

  // This smart pointer will hold the different population types
  boost::shared_ptr<GEvolutionaryAlgorithm> pop_ptr;

  // Create the actual populations
  switch (theAppParams.getParallelizationMode()) {
  case 0: // Serial execution
    {
      // Create an empty population
      pop_ptr = boost::shared_ptr<GEvolutionaryAlgorithm>(new GEvolutionaryAlgorithm());
    }
    break;

  case 1: // Multi-threaded execution
    {
      // Create the multi-threaded population
      boost::shared_ptr<GMultiThreadedEA> popPar_ptr(new GMultiThreadedEA());

      // Population-specific settings
      popPar_ptr->setNThreads(theAppParams.getNEvaluationThreads());

      // Assignment to the base pointer
      pop_ptr = popPar_ptr;
    }
    break;

  case 2: // Networked execution (server-side)
    {
      // Create a network consumer and enrol it with the broker
      boost::shared_ptr<gc::GAsioTCPConsumerT<gg::GIndividual> > gatc(new gc::GAsioTCPConsumerT<gg::GIndividual>(theAppParams.getPort()));
      gatc->setSerializationMode(theAppParams.getSerMode());
      GINDIVIDUALBROKER->enrol(gatc);

      // Create the actual broker population
      boost::shared_ptr<GBrokerEA> popBroker_ptr(new GBrokerEA());
      popBroker_ptr->setWaitFactor(theAppParams.getWaitFactor());

      // Assignment to the base pointer
      pop_ptr = popBroker_ptr;
    }
    break;
  }

  //----------------------------------------------------------------------------------------------
  // Now we have suitable populations and can fill them with data

  // Add individuals to the population
  for(std::size_t p = 0 ; p<theAppParams.getNParents(); p++) {
    pop_ptr->push_back(parentIndividuals[p]);
  }

  // Specify some general population settings
  pop_ptr->setDefaultPopulationSize(theAppParams.getPopulationSize(),theAppParams.getNParents());
  pop_ptr->setMaxIteration(theAppParams.getMaxIterations());
  pop_ptr->setMaxTime(boost::posix_time::minutes(theAppParams.getMaxMinutes()));
  pop_ptr->setReportIteration(theAppParams.getReportIteration());
  pop_ptr->setRecombinationMethod(theAppParams.getRScheme());
  pop_ptr->setSortingScheme(theAppParams.getSmode());
  
  // Do the actual optimization
  pop_ptr->optimize();

  for(std::size_t i=0; i<pop_ptr->size(); i++) {
    std::cout << i << ": " << pop_ptr->at(i)->fitness() << std::endl;
  }

  //----------------------------------------------------------------------------------------------

  boost::shared_ptr<gp::GOmegaPiIndividual> bestIndividual_ptr=pop_ptr->getBestIndividual<gp::GOmegaPiIndividual>();
  assert(bestIndividual_ptr->getFitParams(finalFitParm));
  finalOmegaPiLh=bestIndividual_ptr->omegaPiLhPtr();

  Info << "GenEvA done.\n" << endmsg;

  return bl::tribool(true);
}

/************************************************************************************************/
/**
 * Starts the actual optimization process, using the Minuit library
 *
 * @return A boost::logic::tribool indicating whether the optimization has succeeded
 */
bl::tribool Minuit(
     ApplicationParameter &theAppParams
     , OmegaPiData::fitParamVal &finalFitParm
     , boost::shared_ptr<AbsOmegaPiLh> &finalOmegaPiLh
) {
  Info << "Minuit fit starts.\n" << endmsg;  

  // get pbarpToOmegaPi0States pointer back

  boost::shared_ptr<const pbarpToOmegaPi0States> theOmegaPi0StatesPtr=finalOmegaPiLh->omegaPi0States();

  theOmegaPi0StatesPtr->print(std::cout);

  rm::MOmegaPiFcn mOmegaPiFcn(finalOmegaPiLh);
  rm::MnUserParameters upar;
  minuitStartParam theUserPar;

  if(theAppParams.getAppExecMode() == ApplicationParameter::Minuit)
    {
      if (!theAppParams.getPathStartParamFile().empty())
	{
	  ifstream theFile(theAppParams.getPathStartParamFile().c_str());   
	  if (theFile)
	    {
	      Info << "Using start parameters from file " << theAppParams.getPathStartParamFile() << "\n" << endmsg;
	      theUserPar.ParseStream(theFile);
	      finalOmegaPiLh->setMnUsrParams(upar,theUserPar);
	    }
	  else 
	    {
	      Info << "Start parameter file " << theAppParams.getPathStartParamFile() << " not found!" << endmsg;
	      return false;
	    }
	}
      else finalOmegaPiLh->setMnUsrParams(upar);
    }
  else if(theAppParams.getAppExecMode() == ApplicationParameter::GenToMinuit) finalOmegaPiLh->setMnUsrParams(upar,finalFitParm);
   
  rm::MnMigrad migrad(mOmegaPiFcn, upar);
  Info <<"start migrad "<< endmsg;
  rm::FunctionMinimum min = migrad();

  if(!min.IsValid()) {
    //try with higher strategy
    Info <<"FM is invalid, try with strategy = 2."<< endmsg;
    rm::MnMigrad migrad2(mOmegaPiFcn, min.UserState(), rm::MnStrategy(2));
    min = migrad2();
  }

  rm::MnUserParameters finalUsrParameters=min.UserParameters();
  const std::vector<double> finalParamVec=finalUsrParameters.Params();

  cout << endl << "Errors:" << endl;
  const vector<double> finalParamErrors=finalUsrParameters.Errors();
  
  vector<double>::const_iterator it;
  int i;
  for (it = finalParamErrors.begin(), i=0; it < finalParamErrors.end(); it++,i++) 
    {
      std::string strName =finalUsrParameters.GetName(i);
      cout << "Name:" << strName << " ";
      cout << "Index:" << finalUsrParameters.Index(strName) << " ";
      cout << *it << endl;
    }
  
  finalOmegaPiLh->getFitParamVal(finalFitParm, finalParamVec); 
  
  Info << "Minuit done.\n" << endmsg;

  return bl::tribool(true);
}

/************************************************************************************************/
/**
 *
 */
bool QAmode(boost::shared_ptr<const OmegaPiEventList> &theOmegaPiEventPtr,
		   boost::shared_ptr<const pbarpToOmegaPi0States> &pbarpToOmegaPi0StatesPtr,
		   ApplicationParameter &theAppParams,
		   OmegaPiData::fitParamVal &finalFitParm,
		   boost::shared_ptr<AbsOmegaPiLh> &finalOmegaPiLh
		   )
{
  Info << "QA mode start.\n" << endmsg;
    

  // get pbarpToOmegaPi0States pointer back
  boost::shared_ptr<const pbarpToOmegaPi0States> theOmegaPi0StatesPtr=finalOmegaPiLh->omegaPi0States();  
    
  theOmegaPi0StatesPtr->print(std::cout);


  rm::MOmegaPiFcn mOmegaPiFcn(finalOmegaPiLh);

  rm::MnUserParameters upar;
  minuitStartParam theUserPar;
    

  if (!theAppParams.getPathStartParamFile().empty())
    {
      ifstream theFile(theAppParams.getPathStartParamFile().c_str());   
      if (theFile)
	{
	  Info << "Using start parameters from file " << theAppParams.getPathStartParamFile() << "\n" << endmsg;
	  theUserPar.ParseStream(theFile);
	  finalOmegaPiLh->setMnUsrParams(upar,theUserPar);
	}
      else 
	{
	  Info << "Start parameter file " << theAppParams.getPathStartParamFile() << " not found!" << endmsg;
	  return false;
	}
    }
  else 
    {
      Info << "Start parameter file has to be set!" << endmsg;
      return false;
    }
      
//   mOmegaPiFcn.setFitParamVal(finalFitParm, upar.Params());
  finalOmegaPiLh->getFitParamVal(finalFitParm, upar.Params());	
  
  Info << "QA mode done.\n" << endmsg;
	
  return true;
}

/************************************************************************************************/
/**
 *
 */
bool calcSpinDensity(ApplicationParameter &theAppParams, 
		     boost::shared_ptr<AbsOmegaPiLh> finalOmegaPiLh,
		     OmegaPiData::fitParamVal &finalFitParm)
{
  Info << "Starting spin density calculation." << "\n" << endmsg;
  OmegaPiData::fitParamVal theParamVal;
        
  rm::MOmegaPiFcn mOmegaPiFcn(finalOmegaPiLh);
  rm::MnUserParameters upar;
  minuitStartParam theUserPar;

  if (!theAppParams.getPathStartParamFile().empty())
    {
      ifstream theFile(theAppParams.getPathStartParamFile().c_str());   
      if (theFile)
	{
	  Info << "Using start parameters from file " << theAppParams.getPathStartParamFile() << "\n" << endmsg;
	  theUserPar.ParseStream(theFile);
	  finalOmegaPiLh->setMnUsrParams(upar,theUserPar);
	}
      else 
	{
	  Info << "Start parameter file " << theAppParams.getPathStartParamFile() << " not found!" << endmsg;
	  return false;
	}
    }
  else 
    {
      Info << "Start parameter file has to be set!" << endmsg;
      return false;
    }

  finalOmegaPiLh->getFitParamVal(theParamVal, upar.Params());       
  
  Info << "Using following fit parameter:\n" << endmsg;
  finalOmegaPiLh->printFitParams(std::cout, theParamVal);

  std::ostringstream theSpinDensityRootFile;
        
  if (theAppParams.getCalcAllSpindensity())
    {
      Info << "Calculating all spin density elements.\n" << endmsg;
      theSpinDensityRootFile << "./" << theAppParams.getName() << "SpinDensity" << "_jmax" << theAppParams.getJMax() << "_mom" << theAppParams.getPbarMom() << ".root";

      SpinDensityHist theSpinDensityHist(theSpinDensityRootFile.str(), finalOmegaPiLh, theParamVal);
      theSpinDensityHist.createHistograms();
    }
  else
    {
      Info << "Calculating spin density elements for M=" << theAppParams.getM() << " M'=" << theAppParams.getM_() << ".\n" << endmsg;
      theSpinDensityRootFile << "./" << theAppParams.getName() << "SpinDensity_M1" << theAppParams.getM() << "_M2" << theAppParams.getM_() << "_jmax" << theAppParams.getJMax() << "_mom" << theAppParams.getPbarMom() << ".root";
       
      SpinDensityHist theSpinDensityHist(theSpinDensityRootFile.str(), finalOmegaPiLh, theParamVal);
      theSpinDensityHist.createHistogram(theAppParams.getM(),theAppParams.getM_());
    }
  Info << "Spin density calculation done." << "\n" << endmsg;
  return true;
}

/************************************************************************************************/
/**
 *
 */
void removeEvents(EventList &piOmegaEventsData, int nEventsToRemove, bool bRemoveFromEnd)
{
  if (nEventsToRemove != 0)
    {
      std::stringstream strMsg;
      strMsg << "Removing " << nEventsToRemove << " events";
      if (bRemoveFromEnd)
	{
	  strMsg << " from the end of the event list."  << endl;
	  piOmegaEventsData.removeEvents(piOmegaEventsData.size()-nEventsToRemove,piOmegaEventsData.size());
	}
      else
	{
	  strMsg << " from the beginning of the event list." << endl;
	  piOmegaEventsData.removeEvents(0,piOmegaEventsData.size()-nEventsToRemove);
	}
      Info << strMsg.str() << endmsg;
      Info << "\nEvent list now contains " << piOmegaEventsData.size() << " events. Each event has "
	   <<  piOmegaEventsData.nextEvent()->size() << " final state particles.\n" << endmsg;
    }
}

/************************************************************************************************/
/**
 * Sets the error logging mode
 * 
 * @param erlMode The desired error logging mode
 */
void setErrLogMode( const ApplicationParameter::enErrLogMode& erlMode ) {
  switch(erlMode) {
  case ApplicationParameter::debug :
    ErrLogger::instance()->setLevel(log4cpp::Priority::DEBUG);
    break;
  case ApplicationParameter::trace :
    ErrLogger::instance()->setLevel(log4cpp::Priority::INFO);
    break;
  case ApplicationParameter::routine :
    ErrLogger::instance()->setLevel(log4cpp::Priority::INFO);
    break;
  case ApplicationParameter::warning :
    ErrLogger::instance()->setLevel(log4cpp::Priority::WARN);
    break;
  case ApplicationParameter::error :
    ErrLogger::instance()->setLevel(log4cpp::Priority::ERROR);
    break;
  case ApplicationParameter::alert :
    ErrLogger::instance()->setLevel(log4cpp::Priority::ALERT);
    break;
  default: 
    ErrLogger::instance()->setLevel(log4cpp::Priority::DEBUG);
  }
}

/************************************************************************************************/
/**
 * Inform the audience about the execution mode
 *
 * @param execMode The chosen execution mode
 */
void emitExecutionMode( const ApplicationParameter::enExecMode& execMode ) {
  switch(execMode) {
    case ApplicationParameter::GenEvA:
      Info << "Using minimization algorithm: " << "GenEvA" << "\n" << endmsg;
      break;
    case ApplicationParameter::Minuit:
      Info << "Using minimization algorithm: " << "Minuit" << "\n" << endmsg;
      break;
    case ApplicationParameter::GenToMinuit:
      Info << "Using minimization algorithm: " << "First GenEvA then Minuit with final fit parameters from GenEvA" << "\n" << endmsg;
      break;
    case ApplicationParameter::SpinDensity:
      Info << "Calculating Spin Density" << "\n" << endmsg;
      break;
    case ApplicationParameter::QAmode:
      Info << "using QA mode" << "\n" << endmsg;
      break;
  }
}

/************************************************************************************************/
/**
 * The main function.
 */
int main(int argc, char **argv)
{
  cout << "Combined GenEvA and Minuit fit program " << argv[0] << endl << endl;

  // Parse the command line
  static ApplicationParameter theAppParams(argc, argv);
    
  // Random numbers are our most valuable good. Set the number of threads
  GRANDOMFACTORY->setNProducerThreads(theAppParams.getNProducerThreads());
  GRANDOMFACTORY->setArraySize(theAppParams.getArraySize());

  // Set the desired error logging mode
  setErrLogMode(theAppParams.getErrLogMode());

  // Inform the audience about the execution mode
  emitExecutionMode(theAppParams.getAppExecMode());

  std::string theCfgFile = theAppParams.getConfigFile();
  Info << "The path to config file is " << theCfgFile << "\n" << endmsg;
  Info << "Maximum spin content: " << theAppParams.getJMax() << endmsg;
  Info << "pbar momentum: " << theAppParams.getPbarMom()   << endmsg;
      
  std::string piomegaDatFile;
  std::string piomegaMcFile; 
  
  constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPi/data/510_",theAppParams.getPbarMom(),piomegaDatFile);
  if (checkFileExist(piomegaDatFile))
    {
      Info << "Using Data file " << piomegaDatFile << endmsg;
    }
  else
    {
      Alert <<"Data file for pbarMom= " << theAppParams.getPbarMom() << " not available !"  << endmsg;
      Alert << "File " << piomegaDatFile << " is missing !" << endmsg;
      exit(1);
    }

  constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPi/data/mc510_",theAppParams.getPbarMom(),piomegaMcFile);
  if (checkFileExist(piomegaMcFile))
    {
      Info << "Using Monte Carlo file " << piomegaMcFile << endmsg;
    }
  else
    {
      Alert <<"Monte Carlo file for pbarMom= " << theAppParams.getPbarMom() << "not available !" << endmsg;
      Alert << "File " << piomegaMcFile << " is missing !" << endmsg;
      exit(1);
    }
  
  Info << "data file: " << piomegaDatFile << endmsg;
  Info << "mc file: " << piomegaMcFile << endmsg;

  ParticleTable pTable;
  PdtParser parser;
  std::string pdtFile(theAppParams.getSourcePath()+"/Particle/pdt.table");
  if (!parser.parse(pdtFile, pTable)) {
    Alert << "Error: could not parse " << pdtFile << endmsg;
    exit(1);
  }

  std::vector<std::string> fileNames;
  fileNames.push_back(piomegaDatFile);
  CBElsaReader eventReader(fileNames, 3, 0); 
  EventList piOmegaEventsData;
  eventReader.fillAll(piOmegaEventsData);

  Info << "\nFile has " << piOmegaEventsData.size() << " events. Each event has "
       <<  piOmegaEventsData.nextEvent()->size() << " final state particles.\n" << endmsg;

  if(theAppParams.getNumEventsRemove() > 0) {
    removeEvents(piOmegaEventsData,theAppParams.getNumEventsRemove(),theAppParams.getRemoveEventsFromEnd());
  }
  else if(theAppParams.getNumEventsRed() > 0) {
    int nEventsRed = theAppParams.getNumEventsRed();
    int nEventListSize = piOmegaEventsData.size();
    
    Info << "Reducing number of the events in the event list to " 
	 << nEventsRed << " events." << endmsg;
    
    int nNewSize = nEventListSize-nEventsRed;
    
    if (nNewSize > 0 && nNewSize != nEventListSize) {
      removeEvents(piOmegaEventsData, nNewSize, theAppParams.getRemoveEventsFromEnd());
    }
    else {
      Alert << "Could not reduce number of the parameters because number to which they should be"
	    << " reduced " << "(" << nEventsRed << ")" 
	    <<" to is to great or same as the actual number of events in the list.\n" 
	    << endmsg;
      exit(-1);
    }
  }
      
  if (!piOmegaEventsData.findParticleTypes(pTable)) {
    Warning << "could not find all particles" << endmsg;
  }
  piOmegaEventsData.rewind();


  Event* anEvent;
  int evtCount = 0;
  while ((anEvent = piOmegaEventsData.nextEvent()) != 0 && evtCount < 20) {
    Info << "\n" 
	 << *(anEvent->p4(0)) << "\tm = " << anEvent->p4(0)->Mass() << "\n"
	 << *(anEvent->p4(1)) << "\tm = " << anEvent->p4(1)->Mass() << "\n"
	 << *(anEvent->p4(2)) << "\tm = " << anEvent->p4(2)->Mass() << "\n"
	 << endmsg;
    ++evtCount;
  }
  piOmegaEventsData.rewind();

  std::vector<std::string> fileNamesMc;

  fileNamesMc.push_back(piomegaMcFile);
  CBElsaReader eventReaderMc(fileNamesMc, 3, 0); 
  EventList piOmegaEventsMc;
  eventReaderMc.fillAll(piOmegaEventsMc);
  piOmegaEventsMc.rewind();

  boost::shared_ptr<const OmegaPiEventList> theOmegaPiEventPtr(new OmegaPiEventList(piOmegaEventsData, piOmegaEventsMc, theAppParams.getJMax(),  theAppParams.getPbarMom()));
  boost::shared_ptr<pbarpStates> pbarpStatesPtr(new pbarpStates(theAppParams.getJMax()));
  boost::shared_ptr<const pbarpToOmegaPi0States> pbarpToOmegaPi0StatesPtr(new pbarpToOmegaPi0States(pbarpStatesPtr));
  OmegaPiData::fitParamVal finalFitParm;


  boost::shared_ptr<AbsOmegaPiLh> finalOmegaPiLh;
  if (theAppParams.getLhMode()=="OmegaPiLhGamma"){
    finalOmegaPiLh = boost::shared_ptr<AbsOmegaPiLh>(new OmegaPiLhGamma(theOmegaPiEventPtr, pbarpToOmegaPi0StatesPtr));  
  }
  else if (theAppParams.getLhMode()=="OmegaPiLhGammaBw"){
    finalOmegaPiLh = boost::shared_ptr<AbsOmegaPiLh>(new OmegaPiLhGammaBw(theOmegaPiEventPtr, pbarpToOmegaPi0StatesPtr));  
  }
  else{
    Alert <<"LhMode " << theAppParams.getLhMode() << " doesn't exist !!! " << endmsg;
    exit(1);
  }

  Info << "\n using Lh mode:"; 
  finalOmegaPiLh->print(std::cout);

  //---------------------------------------------------------------------------
  // The actual optimization
  bl::tribool bExecFinish=bl::tribool(false);
  switch(theAppParams.getAppExecMode()) {
  case ApplicationParameter::GenEvA:
    {
      bExecFinish = GenEvA(theOmegaPiEventPtr,pbarpToOmegaPi0StatesPtr,theAppParams,finalFitParm,finalOmegaPiLh);
      if(bl::indeterminate(bExecFinish)) return 0; // Simply terminate -- this is a client in networked mode that has done its job

      std::ostringstream theParamFilePath;
      theParamFilePath << "./" << theAppParams.getName() << "LastFitParamOmegaPi0Fit_jmax" << theAppParams.getJMax() << "_mom" << theAppParams.getPbarMom() << ".txt";
      ofstream outfile (theParamFilePath.str().c_str());
      //       outfile << PrintFinalFitParam(finalFitParm);
      finalOmegaPiLh->dumpCurrentResult(outfile, finalFitParm);
      outfile.close();
    }
    break;

  case ApplicationParameter::Minuit:
    {
      bExecFinish = Minuit(theAppParams,finalFitParm,finalOmegaPiLh);
      std::ostringstream theParamFilePathMin;
      theParamFilePathMin << "./" << theAppParams.getName() << "LastFitParamOmegaPi0Fit_jmax" << theAppParams.getJMax() << "_mom" << theAppParams.getPbarMom() << ".txt";
      ofstream outfileMin (theParamFilePathMin.str().c_str());
    }
    break;

  case ApplicationParameter::GenToMinuit:
    {
      bExecFinish = GenEvA(theOmegaPiEventPtr,pbarpToOmegaPi0StatesPtr,theAppParams,finalFitParm,finalOmegaPiLh);
      Info << "Obtained NLL with GenEvA:\n" << finalOmegaPiLh->calcLogLh(finalFitParm) << endmsg;
      Info << "\n and fit parameters are:\n" << endmsg;
      finalOmegaPiLh->printFitParams(std::cout, finalFitParm);
      std::ostringstream theParamFilePathGen;
      theParamFilePathGen << "./" << theAppParams.getName() << "GenevaLastFitParamOmegaPi0Fit_jmax" << theAppParams.getJMax() << "_mom" << theAppParams.getPbarMom() << ".txt";
      ofstream outfileGen (theParamFilePathGen.str().c_str());
      finalOmegaPiLh->dumpCurrentResult(outfileGen, finalFitParm);
      outfileGen.close();
      
      if(bl::indeterminate(bExecFinish)) return 0; // Simply terminate -- this is a client in networked mode that has done its job
      bExecFinish = Minuit(theAppParams,finalFitParm,finalOmegaPiLh);
      std::ostringstream theParamFilePathMin;
      theParamFilePathMin << "./" << theAppParams.getName() << "LastFitParamOmegaPi0Fit_jmax" << theAppParams.getJMax() << "_mom" << theAppParams.getPbarMom() << ".txt";
      ofstream outfileMin (theParamFilePathMin.str().c_str());
      finalOmegaPiLh->dumpCurrentResult(outfileMin, finalFitParm);
      outfileMin.close();
    }
    break;

  case ApplicationParameter::SpinDensity:
    {
      if (!calcSpinDensity(theAppParams, finalOmegaPiLh, finalFitParm)) exit(1);
    }
    break;

  case ApplicationParameter::QAmode:
    {
      bExecFinish = QAmode(theOmegaPiEventPtr,pbarpToOmegaPi0StatesPtr,theAppParams,finalFitParm,finalOmegaPiLh);
    }
    break;

  default:
    {
      cerr << "Error unknown execution mode selected!" << endl;
    }
  }

  //---------------------------------------------------------------------------
  // Let the audience know
  if (bExecFinish) {
      cout << endl;
      Info << "Fit results:\n" << endmsg;
      Info << "Final logLH=" << finalOmegaPiLh->calcLogLh(finalFitParm) << "\n" << endmsg;
      Info << "Final fit parameters:\n" << endmsg;
//       printFitParameters(pbarpToOmegaPi0StatesPtr, finalFitParm);
      finalOmegaPiLh->printFitParams(std::cout, finalFitParm);
      std::ostringstream theRootFilePath;
      theRootFilePath << "./" << theAppParams.getName() << "OmegaPi0Fit_jmax" << theAppParams.getJMax() << "_mom" << theAppParams.getPbarMom() << ".root";
      OmegaPiHist theHistogrammer(finalOmegaPiLh,finalFitParm,theRootFilePath.str());
  }

  //---------------------------------------------------------------------------

  return 0;
}

/************************************************************************************************/