#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <cstdlib>
#include <algorithm>

#include "TROOT.h"
#include "TFile.h"
#include "TTree.h"
#include "TBits.h"
#include "TObjString.h"
#include "TGeoManager.h"

using namespace std;

namespace {
  vector<int> inputNeutrinoTypes; // vector of PDGs of allowed neutrinos
  vector<pair<double, double> > inputNeutrinoEnergyRanges; // vector of allowed (min, max) neutrino energies
  vector<int> inputNuclearTargets; // vector of allowed nuclear targets
  vector<string> geoTargets; // geometric target pathnames
  vector<int> reactionCodes; // vector of allowed reactions
  map<int, int> outputParticles; // map of pdg:number of output particles
  int countLeptons = -1;

  vector<string> inputFileNames;
  string outputFileName = "";
  string treeName = "";

  bool hasSelectionCriteria = false;
  bool invert = false;

  const unsigned int kNPmax = 100;

  bool verbose = false;

  //-- define the rootracker tree branches
  TBits*      brEvtFlags = 0;             // generator-specific event flags
  TObjString* brEvtCode = 0;              // generator-specific string with 'event code'
  int         brEvtNum;                   // event num.
  double      brEvtXSec;                  // cross section for selected event (1E-38 cm2)
  double      brEvtDXSec;                 // cross section for selected event kinematics (1E-38 cm2 /{K^n})
  double      brEvtWght;                  // weight for that event
  double      brEvtProb;                  // probability for that event (given cross section, path lengths, etc)
  double      brEvtVtx[4];                // event vertex position in detector coord syst (SI)
  int         brStdHepN;                  // number of particles in particle array 
  // > stdhep-like particle array:
  int         brStdHepPdg   [kNPmax];     // pdg codes (& generator specific codes for pseudoparticles)
  int         brStdHepStatus[kNPmax];     // generator-specific status code
  int         brStdHepRescat[kNPmax];     // Hadron transport model - specific rescattering code
  double      brStdHepX4    [kNPmax][4];  // 4-x (x, y, z, t) of particle in hit nucleus frame (fm)
  double      brStdHepP4    [kNPmax][4];  // 4-p (px,py,pz,E) of particle in LAB frame (GeV)
  double      brStdHepPolz  [kNPmax][3];  // polarization vector
  int         brStdHepFd    [kNPmax];     // first daughter
  int         brStdHepLd    [kNPmax];     // last  daughter 
  int         brStdHepFm    [kNPmax];     // first mother
  int         brStdHepLm    [kNPmax];     // last  mother
  // > neutrino parent info (passed-through from the beam-line MC / quantities in 'jnubeam' units)
  int         brNuParentPdg;              // parent hadron pdg code
  int         brNuParentDecMode;          // parent hadron decay mode
  double      brNuParentDecP4 [4];        // parent hadron 4-momentum at decay 
  double      brNuParentDecX4 [4];        // parent hadron 4-position at decay
  double      brNuParentProP4 [4];        // parent hadron 4-momentum at production
  double      brNuParentProX4 [4];        // parent hadron 4-position at production
  int         brNuParentProNVtx;          // parent hadron vtx id
  // variables added since 10a flux compatibility changes
  int         brNuFluxEntry;              // Entry number from flux file
  int         brNuIdfd;                   // Detector location id
  double      brNuCospibm;                // Cosine of the angle between the parent particle direction and the beam direction
  double      brNuCospi0bm;               // Same as above except at the production of the parent particle
  int         brNuGipart;                 // Primary particle ID
  float       brNuGpos0[3];               // Primary particle starting point
  float       brNuGvec0[3];               // Primary particle direction at the starting point
  float       brNuGamom0;                 // Momentum of the primary particle at the starting point
  // > etc 
  int         brNeutCode;                 // > NEUT-like reaction code for the GENIE event

  TTree*      inputTree;
  TTree*      outputTree;

  TFile*      inputFile;
  TFile*      outputFile;

  TGeoManager* geoManager = 0;
}

bool ParseArgs(int argc, char** argv);
bool ParseInputNeutrinos(const char* optstr);
bool ParseInputNeutrinoEnergies(const char* optstr);
bool ParseInputNuclearTargets(const char* optstr);
bool ParseReactionCode(const char* optstr);
void ParseGeoTargets(const char* optstr);
bool ParseOutputParticles(const char* optstr);

int MatchParticlePdg(const string& particle);
int AtomicNumber(const string& atom);

inline bool SelectThisEvent();
bool SelectInput();
bool SelectInteraction();
bool CountParticles();

void PrintUsage();
void PrintHelp();
void PrintOptions();

int main(int argc, char** argv) {
  if(!ParseArgs(argc, argv)) {
    cout << "Invalid option(s), exiting" << endl;
    return 1;
  }
  while(argc - optind > 0) {
    const char* filename = argv[optind++];
    TFile tmpfile(filename, "READ");
    if(tmpfile.IsZombie()) {
      cerr << "Error: file " << filename << " does not exist" << endl;
      return 1;
    }
    treeName = "gRooTracker";
    inputTree = dynamic_cast<TTree*>(tmpfile.Get(treeName.c_str()));
    if(!inputTree) {
      treeName = "nRooTracker";
      inputTree = dynamic_cast<TTree*>(tmpfile.Get(treeName.c_str()));
    }
    if(!inputTree) {
      cerr << "Error: RooTracker tree not found in " << filename << endl;
      return 1;
    }
    inputFileNames.push_back(filename);
  }
  if(!(inputFileNames.size() > 0 && outputFileName != "" && hasSelectionCriteria)) {
    PrintUsage();
    return 1;
  }
  
  outputFile = new TFile(outputFileName.c_str(), "RECREATE");
  //outputTree = new TTree(treeName.c_str(),"Event tree rootracker format");
  //-- create the output ROOT tree branches
  /*
  outputTree->Branch("EvtFlags", "TBits",      &brEvtFlags, 32000, 1);           
  outputTree->Branch("EvtCode",  "TObjString", &brEvtCode,  32000, 1);            
  outputTree->Branch("EvtNum",          &brEvtNum,          "EvtNum/I");             
  outputTree->Branch("EvtXSec",         &brEvtXSec,         "EvtXSec/D");            
  outputTree->Branch("EvtDXSec",        &brEvtDXSec,        "EvtDXSec/D");           
  outputTree->Branch("EvtWght",         &brEvtWght,         "EvtWght/D");            
  outputTree->Branch("EvtProb",         &brEvtProb,         "EvtProb/D");            
  outputTree->Branch("EvtVtx",           brEvtVtx,          "EvtVtx[4]/D");             
  outputTree->Branch("StdHepN",         &brStdHepN,         "StdHepN/I");              
  outputTree->Branch("StdHepPdg",        brStdHepPdg,       "StdHepPdg[StdHepN]/I");  
  outputTree->Branch("StdHepStatus",     brStdHepStatus,    "StdHepStatus[StdHepN]/I"); 
  outputTree->Branch("StdHepX4",         brStdHepX4,        "StdHepX4[StdHepN][4]/D"); 
  outputTree->Branch("StdHepP4",         brStdHepP4,        "StdHepP4[StdHepN][4]/D"); 
  outputTree->Branch("StdHepPolz",       brStdHepPolz,      "StdHepPolz[StdHepN][3]/D"); 
  outputTree->Branch("StdHepFd",         brStdHepFd,        "StdHepFd[StdHepN]/I"); 
  outputTree->Branch("StdHepLd",         brStdHepLd,        "StdHepLd[StdHepN]/I"); 
  outputTree->Branch("StdHepFm",         brStdHepFm,        "StdHepFm[StdHepN]/I"); 
  outputTree->Branch("StdHepLm",         brStdHepLm,        "StdHepLm[StdHepN]/I");
  outputTree->Branch("NuParentPdg",     &brNuParentPdg,     "NuParentPdg/I");       
  outputTree->Branch("NuParentDecMode", &brNuParentDecMode, "NuParentDecMode/I");   
  outputTree->Branch("NuParentDecP4",    brNuParentDecP4,   "NuParentDecP4[4]/D");     
  outputTree->Branch("NuParentDecX4",    brNuParentDecX4,   "NuParentDecX4[4]/D");     
  outputTree->Branch("NuParentProP4",    brNuParentProP4,   "NuParentProP4[4]/D");     
  outputTree->Branch("NuParentProX4",    brNuParentProX4,   "NuParentProX4[4]/D");     
  outputTree->Branch("NuParentProNVtx", &brNuParentProNVtx, "NuParentProNVtx/I");   
  outputTree->Branch("NuGipart",        &brNuGipart,        "NuGipart/I");
  outputTree->Branch("NuGpos0",          brNuGpos0,         "NuGpos0[3]/D");
  outputTree->Branch("NuGvec0",          brNuGvec0,         "NuGvec0[3]/D");
  outputTree->Branch("NuGamom0",        &brNuGamom0,        "NuGamom0/D");
  // the following branches only exist in GENIE rootracker format
  if(treeName == "gRooTracker") {
    outputTree->Branch("StdHepRescat",     brStdHepRescat,    "StdHepRescat[StdHepN]/I");
    outputTree->Branch("NuFluxEntry",     &brNuFluxEntry,     "NuFluxEntry/I");
    outputTree->Branch("NuIdfd",          &brNuIdfd,          "NuIdfd/I");
    outputTree->Branch("NuCospibm",       &brNuCospibm,       "NuCospibm/D");
    outputTree->Branch("NuCospi0bm",      &brNuCospi0bm,      "NuCospi0bm/D");
    outputTree->Branch("G2NeutEvtCode",   &brNeutCode,        "G2NeutEvtCode/I");
  }
  */

  PrintOptions();
  
  double sumWeight = 0.;

  for(vector<string>::iterator filename = inputFileNames.begin(); filename != inputFileNames.end(); ++filename) {
    cout << "Opening " << *filename << endl;
    inputFile = new TFile(filename->c_str(), "READ");
    inputTree = dynamic_cast<TTree*>(inputFile->Get(treeName.c_str()));

    outputFile->cd();
    outputTree = (TTree*) inputTree->CloneTree(0);

    inputTree->SetBranchAddress("EvtFlags",       &brEvtFlags);
    inputTree->SetBranchAddress("EvtCode",        &brEvtCode);
    inputTree->SetBranchAddress("EvtNum",         &brEvtNum);
    inputTree->SetBranchAddress("EvtXSec",        &brEvtXSec);
    inputTree->SetBranchAddress("EvtDXSec",       &brEvtDXSec);
    inputTree->SetBranchAddress("EvtWght",        &brEvtWght);
    inputTree->SetBranchAddress("EvtProb",        &brEvtProb);
    inputTree->SetBranchAddress("EvtVtx",          brEvtVtx);
    inputTree->SetBranchAddress("StdHepN",        &brStdHepN);
    inputTree->SetBranchAddress("StdHepPdg",       brStdHepPdg);
    inputTree->SetBranchAddress("StdHepStatus",    brStdHepStatus);
    inputTree->SetBranchAddress("StdHepX4",        brStdHepX4);
    inputTree->SetBranchAddress("StdHepP4",        brStdHepP4);
    inputTree->SetBranchAddress("StdHepPolz",      brStdHepPolz);
    inputTree->SetBranchAddress("StdHepFd",        brStdHepFd);
    inputTree->SetBranchAddress("StdHepLd",        brStdHepLd);
    inputTree->SetBranchAddress("StdHepFm",        brStdHepFm);
    inputTree->SetBranchAddress("StdHepLm",        brStdHepLm);
    inputTree->SetBranchAddress("NuParentPdg",    &brNuParentPdg);
    inputTree->SetBranchAddress("NuParentDecMode",&brNuParentDecMode);
    inputTree->SetBranchAddress("NuParentDecP4",   brNuParentDecP4);
    inputTree->SetBranchAddress("NuParentDecX4",   brNuParentDecX4);
    inputTree->SetBranchAddress("NuParentProP4",   brNuParentProP4);
    inputTree->SetBranchAddress("NuParentProX4",   brNuParentProX4);
    inputTree->SetBranchAddress("NuParentProNVtx",&brNuParentProNVtx);   
    inputTree->SetBranchAddress("NuGipart",       &brNuGipart);
    inputTree->SetBranchAddress("NuGpos0",         brNuGpos0);
    inputTree->SetBranchAddress("NuGvec0",         brNuGvec0);
    inputTree->SetBranchAddress("NuGamom0",       &brNuGamom0);
    // the following branches only exist in GENIE rootracker format
    if(treeName == "gRooTracker") {
      inputTree->SetBranchAddress("StdHepRescat",    brStdHepRescat);
      inputTree->SetBranchAddress("NuFluxEntry",    &brNuFluxEntry);
      inputTree->SetBranchAddress("NuIdfd",         &brNuIdfd);
      inputTree->SetBranchAddress("NuCospibm",      &brNuCospibm);
      inputTree->SetBranchAddress("NuCospi0bm",     &brNuCospi0bm);
      inputTree->SetBranchAddress("G2NeutEvtCode",  &brNeutCode);
    }




    int nInput = 0, nOutput = 0;

    Long64_t N = inputTree->GetEntries();
    for(Long64_t i = 0; i < N; ++i) {
      inputTree->GetEntry(i);
      nInput ++;
      if(SelectThisEvent()) {
        outputTree->Fill();
        nOutput++;
      }
      if(verbose && nInput % 5000 == 0) {
        cout << nInput << " events processed, " << nOutput << " events selected." << std::endl;
      }
    }
    
    sumWeight += inputTree->GetWeight();
    
    cout << nInput << " events processed, " << nOutput << " events selected." << std::endl;
    
    inputFile->Close();
    delete inputFile;
    
    

  }
  outputTree->SetWeight(sumWeight);

  outputFile->cd();
  outputTree->Write();
  outputFile->Write();
  outputFile->Close();

  return 0;
}

void PrintUsage() {
  cout << "Usage: " << endl;
  cout << "RunOACherryPicker.exe [options] input.root -o output.root" << endl;
}

void PrintHelp() {
  cout << "Options: " << endl;
  cout << "   -h:" << endl
       << "        show this help message" << endl;
  cout << "   -n [12,14,16,-12,-14,-16]:" << endl
       << "        select by input neutrino (pdg code)" << endl;
  cout << "   -e [elow,ehigh]:" << endl
       << "        select by neutrino energy range [elow,ehigh]" << endl;
  cout << "   -t [aaa]:" << endl
       << "        select by target nucleus atomic number or name (e.g. \"6\" or \"c\" or \"carbon\")" << endl;
  cout << "   -r [nnn]:" << endl
       << "        select by NEUT interaction code" << endl;
  cout << "   -g [geometry.root] -l [A/B/C]:" << endl
       << "        select by interaction location, e.g. \"FGD/Scint\" for events in the active volumes of either FGD" << endl;
  cout << "   -c [name:n]:" << endl
       << "        select by outgoing particle (pdg code or name) counts, e.g. \"pi0:1\" for events with 1 outgoing pi0" << endl;
  cout << "   -i:" << endl
       << "        invert the selection criteria" << endl;
  cout << "   -v:" << endl
       << "        verbose mode" << endl;
  cout << "Please see the documentation for more explanation and examples." << endl;
}

bool SelectThisEvent() {
  bool select = SelectInput() && SelectInteraction() && CountParticles();
  if(invert) select = !select;
  return select;
}

bool SelectInput() {
  bool select = true;
  if(inputNeutrinoTypes.size() > 0) {
    for(int iparticle = 0; iparticle < brStdHepN; ++iparticle) {
      if(brStdHepStatus[iparticle] != 0) continue;
      int pdg = brStdHepPdg[iparticle];
      if(pdg == 12 || pdg == 14 || pdg == 16 || pdg == -12 || pdg == -14 || pdg == -16) {
        if(find(inputNeutrinoTypes.begin(), inputNeutrinoTypes.end(), pdg) == inputNeutrinoTypes.end()) {
          select = false;
        }
      }
    }
  }
  if(select && inputNeutrinoEnergyRanges.size() > 0) {
    for(int iparticle = 0; iparticle < brStdHepN; ++iparticle) {
      if(brStdHepStatus[iparticle] != 0) continue;
      int pdg = brStdHepPdg[iparticle];
      double E = brStdHepP4[iparticle][3];
      if(pdg == 12 || pdg == 14 || pdg == 16 || pdg == -12 || pdg == -14 || pdg == -16) {
        select = false;
        for(vector<pair<double, double> >::const_iterator e = inputNeutrinoEnergyRanges.begin(); e != inputNeutrinoEnergyRanges.end(); ++e) {
          if(E >= e->first && E <= e->second) {
            select = true;
          }
        }
      }
    }
  }
  if(select && inputNuclearTargets.size() > 0) {
    for(int iparticle = 0; iparticle < brStdHepN; ++iparticle) {
      if(brStdHepStatus[iparticle] != 0) continue;
      int pdg = brStdHepPdg[iparticle];
      if((pdg > 1000000000 && pdg < 1999999999)) {
        int atom = (pdg / 10000) % 1000;
        if(find(inputNuclearTargets.begin(), inputNuclearTargets.end(), atom) == inputNuclearTargets.end()) {
          select = false;
        }
      }
    }
  }
  return select;
}

bool RecursiveFindInGeoPath(const string& geoPath, const string& pathToFind) {
  //cout << "RecursiveFindInGeoPath: Looking for \"" << pathToFind << "\" in \"" << geoPath << "\"" << endl;
  if(pathToFind.find("/") == string::npos) {
    if(geoPath.find(pathToFind) != string::npos) {
      //cout << "... found!" << endl;
      return true;
    }
    else {
      //cout << "... not found" << endl;
      return false;
    }
  }
  else {
    string preSlash, postSlash, pathPostSlash;
    preSlash = pathToFind.substr(0, pathToFind.find("/"));
    if(geoPath.find(preSlash) == string::npos) {
      //cout << "... not found" << endl;
      return false;
    }
    postSlash = pathToFind.substr(pathToFind.find("/") + 1);
    pathPostSlash = geoPath.substr(geoPath.find("/", geoPath.find(preSlash)) + 1);
    //cout << "... found \"" << preSlash << "\" and now recursing ..." << endl;
    return RecursiveFindInGeoPath(pathPostSlash, postSlash);
  }
}

bool SelectInteraction() {
  bool select = true;
  if(reactionCodes.size() > 0) {
    int reactionCode = brNeutCode;
    if(reactionCode == 0) {
      istringstream rc(brEvtCode->GetString().Data());
      rc >> reactionCode;
    }
    select = false;
    for(vector<int>::const_iterator c = reactionCodes.begin(); c != reactionCodes.end(); ++c) {
      if(reactionCode == *c) {
        select = true;
      }
    }
  }
  if(select && geoManager && geoTargets.size() > 0) {
    double x, y, z;
    x = brEvtVtx[0] * 1000;
    y = brEvtVtx[1] * 1000;
    z = brEvtVtx[2] * 1000;
    geoManager->FindNode(x, y, z);
    string path = geoManager->GetPath();
    select = false;
    for(vector<string>::const_iterator p = geoTargets.begin(); p != geoTargets.end(); ++p) {
      if(RecursiveFindInGeoPath(path, *p)) {
        select = true;
      }
    }
  }
  return select;
}

bool CountParticles() {
  bool select = true;
  if(outputParticles.size() > 0) {
    map<int, int> particleCount;
    int leptonCount = 0;
    for(map<int, int>::const_iterator p = outputParticles.begin();
          p != outputParticles.end(); ++p) {
      particleCount[p->first] = 0;
    }
    for(int iparticle = 0; iparticle < brStdHepN; ++iparticle) {
      if(brStdHepStatus[iparticle] != 1) continue;
      int pdg = brStdHepPdg[iparticle];
      if(particleCount.find(pdg) != particleCount.end()) {
        particleCount[pdg] ++;
      }
      if(pdg == 11 || pdg == 13 || pdg == 15 || pdg == -11 || pdg == -13 || pdg == -15) {
        leptonCount ++;
      }
    }
    for(map<int, int>::const_iterator p = outputParticles.begin();
          p != outputParticles.end(); ++p) {
      if(p->second != particleCount[p->first]) {
        select = false;
        break;
      }
    }
    if(select && countLeptons >= 0 && leptonCount != countLeptons) select = false; 
  }
  return select;
}

bool ParseArgs(int argc, char** argv) {
  for(;;) {
    int c = getopt(argc, argv, "hio:g:c:n:t:l:r:e:v");
    if(c < 0) break;
    switch(c) {
    case 'h':
      PrintUsage();
      PrintHelp();
      exit(0);
      break;
    case 'o':
      {
        outputFileName = optarg;
      }
      break;
    case 'g':
      {
        geoManager = TGeoManager::Import(optarg);
        if(!geoManager) {
          cerr << "Error: No geometry found in " << optarg << endl;
          return false;
        }
      }
      break;
    case 'n':
      if(!ParseInputNeutrinos(optarg)) return false;
      break;
    case 'e':
      if(!ParseInputNeutrinoEnergies(optarg)) return false;
      break;
    case 't':
      if(!ParseInputNuclearTargets(optarg)) return false;
      break;
    case 'r':
      if(!ParseReactionCode(optarg)) return false;
      break;
    case 'l':
      ParseGeoTargets(optarg);
      break;
    case 'c':
      if(!ParseOutputParticles(optarg)) return false;
      break;
    case 'i':
      invert = true;
      break;
    case 'v':
      verbose = true;
      break;
    case '?':
      {
        char x[2] = " ";
        x[0] = optopt;
        std::cout << "Unknown command line option: " << x << std::endl;
      }
      PrintUsage();
      return false;
      break;
    }
  }
  return true;
}

bool ParseInputNeutrinos(const char* optstr) {
  string option = optstr;
  int pdg;
  size_t semicolonpos;
  do {
    pdg = 0;
    semicolonpos = option.find(";");
    istringstream neut(option.substr(0, semicolonpos).c_str());
    neut >> pdg;
    if(pdg == 12 || pdg == 14 || pdg == 16 || pdg == -12 || pdg == -14 || pdg == -16) {
      inputNeutrinoTypes.push_back(pdg);
      hasSelectionCriteria = true;
    }
    else {
      cerr << "Error: Neutrino pdg code \"" << option.substr(0, semicolonpos) << "\" is invalid" << endl;
      return false;
    }
    if(semicolonpos != string::npos) {
      option = option.substr(semicolonpos + 1);
    }
  } while(semicolonpos != string::npos);
  return true;
}

bool ParseInputNeutrinoEnergies(const char* optstr) {
  string option = optstr, energyrange;
  double el, eh;
  size_t semicolonpos, dashpos;
  do {
    semicolonpos = option.find(";");
    energyrange = option.substr(0, semicolonpos);
    dashpos = energyrange.find("-");
    if(dashpos == string::npos || dashpos == 0) {
      cerr << "Error, energy needs to be a range in the format \"Elow-Ehigh\": \"" << energyrange << "\"" << endl;
      return false;
    }
    else {
      el = -1.; eh = -1.;
      istringstream sel(energyrange.substr(0, dashpos).c_str());
      sel >> el;
      istringstream seh(energyrange.substr(dashpos+1).c_str());
      seh >> eh;
      if(el < 0.) {
        cerr << "Low energy must be positive, will use 0." << endl;
        el = 0.;
      }
      if(eh < 0.) {
        cerr << "Error, high energy must be positive" << endl;
        return false;
      }
      else {
        if(eh < el) {
          cerr << "Low energy (" << el << ") must be lower than high energy (" << eh << "), will swap" << endl;
          double tmp = el;
          el = eh;
          eh = tmp;
        }
        inputNeutrinoEnergyRanges.push_back(pair<double, double>(el, eh));
        hasSelectionCriteria = true;
      }
    }
    if(semicolonpos != string::npos) {
      option = option.substr(semicolonpos + 1);
    }
  } while(semicolonpos != string::npos);
  return true;
}

bool ParseInputNuclearTargets(const char* optstr) {
  cout << optstr << endl;
  string option = optstr, atomstring;
  int atom;
  size_t semicolonpos;
  for(unsigned int i = 0; i < option.length(); ++i) {
    option[i] = tolower(option[i]);
  }
  bool exclude = false;
  do {
    semicolonpos = option.find(";");
    atomstring = option.substr(0, semicolonpos);
    if(atomstring.find("exclude") != string::npos) {
      hasSelectionCriteria = true;
      exclude = true;
      for(int i = 1; i <= 109; ++i) {
        inputNuclearTargets.push_back(i);
      }
    }
    else if((atom = AtomicNumber(atomstring)) != 0) {
      if(exclude) {
        vector<int>::iterator t = find(inputNuclearTargets.begin(), inputNuclearTargets.end(), atom);
        if(t != inputNuclearTargets.end()) {
          inputNuclearTargets.erase(t);
        }
      }
      else {
        inputNuclearTargets.push_back(atom);
      }
      hasSelectionCriteria = true;
    }
    else {
      istringstream a(atomstring.c_str());
      a >> atom;
      if(atom != 0) {
        if(exclude) {
          vector<int>::iterator t = find(inputNuclearTargets.begin(), inputNuclearTargets.end(), atom);
          if(t != inputNuclearTargets.end()) {
            inputNuclearTargets.erase(t);
          }
        }
        else {
          inputNuclearTargets.push_back(atom);
        }
        hasSelectionCriteria = true;
      }
      else {
        cerr << "Error, target option \"" << atomstring << "\" is invalid" << endl;
        return false;
      }
    }
    if(semicolonpos != string::npos) {
      option = option.substr(semicolonpos + 1);
    }
  } while(semicolonpos != string::npos);
  return true;
}


bool ParseReactionCode(const char* optstr) {
  string option = optstr;
  size_t semicolonpos;
  int reactionCode;
  do {
    reactionCode = 0;
    semicolonpos = option.find(";");
    istringstream rc(option.substr(0, semicolonpos).c_str());
    rc >> reactionCode;
    if(reactionCode != 0) {
      hasSelectionCriteria = true;
      reactionCodes.push_back(reactionCode);
    }
    else {
      cerr << "Error, invalid reaction code: \"" << option.substr(0, semicolonpos) << "\"" << endl;
      return false;
    }
    if(semicolonpos != string::npos) {
      option = option.substr(semicolonpos + 1);
    }
  } while(semicolonpos != string::npos);
  return true;
}

void ParseGeoTargets(const char* optstr) {
  string option = optstr;
  size_t semicolonpos;
  do {
    semicolonpos = option.find(";");
    geoTargets.push_back(option.substr(0, semicolonpos));
    hasSelectionCriteria = true;
    if(semicolonpos != string::npos) {
      option = option.substr(semicolonpos + 1);
    }
  } while(semicolonpos != string::npos);
}

bool ParseOutputParticles(const char* optstr) {
  string option = optstr, particlestring;
  int number, pdg;
  size_t colonpos, semicolonpos;
  for(unsigned int i = 0; i < option.length(); ++i) {
    option[i] = tolower(option[i]);
  }
  do {
    pdg = 0, number = -1;
    semicolonpos = option.find(";");
    particlestring = option.substr(0, semicolonpos);
    colonpos = particlestring.find(":");
    if(colonpos != string::npos) {
      istringstream n(particlestring.substr(colonpos+1).c_str());
      n >> number;
      if(number >= 0) {
        string particletype = particlestring.substr(0, colonpos);
        if(particletype.find("nucleons") != string::npos) {
          outputParticles[2212] = number;
          outputParticles[2112] = number;
          hasSelectionCriteria = true;
        }
        else if(particletype.find("leptons") != string::npos) {
          countLeptons = number;
          hasSelectionCriteria = true;
        }
        else if(particletype.find("mesons") != string::npos) {
          for(int i = 100; i <= 999; ++i) {
            outputParticles[i] = number;
            outputParticles[-i] = number;
          }
          hasSelectionCriteria = true;
        }
        else if((pdg = MatchParticlePdg(particletype)) != 0) {
          outputParticles[pdg] = number;
            hasSelectionCriteria = true;
        }
        else {
          istringstream p(particletype.c_str());
          p >> pdg;
          if(pdg != 0) {
            outputParticles[pdg] = number;
            hasSelectionCriteria = true;
          }
          else {
            cerr << "Error, particle option \"" << particlestring << "\" is invalid" << endl;
            return false;
          }
        }
      }
      else {
        cerr << "Error, particle option \"" << particlestring << "\" is invalid" << endl;
        return false;
      }
    }
    else {
      cerr << "Error, particle option \"" << particlestring << "\" is invalid" << endl;
      return false;
    }
    if(semicolonpos != string::npos) {
      option = option.substr(semicolonpos + 1);
    }
  } while(semicolonpos != string::npos);
  return true;
}

int MatchParticlePdg(const string& particle) {
  if(particle.find("neutron") != string::npos) {
    return 2112;
  }
  if(particle.find("proton") != string::npos) {
    return 2212;
  }
  if(particle.find("pi0") != string::npos) {
    return 111;
  }
  if(particle.find("pi+") != string::npos || particle.find("pip") != string::npos) {
    return 211;
  }
  if(particle.find("pi-") != string::npos || particle.find("pim") != string::npos) {
    return -211;
  }
  if(particle.find("k0") != string::npos) {
    return 130;
  }
  if(particle.find("k+") != string::npos || particle.find("kp") != string::npos) {
    return 321;
  }
  if(particle.find("k-") != string::npos || particle.find("km") != string::npos) {
    return -321;
  }
  if(particle.find("eta") != string::npos) {
    return 221;
  }
  return 0;
}

int AtomicNumber(const string& atom) {
  if(atom == "h"  || atom == "hydrogen")      return 1; 
  if(atom == "he" || atom == "helium")        return 2; 
  if(atom == "li" || atom == "lithium")       return 3; 
  if(atom == "be" || atom == "beryllium")     return 4; 
  if(atom == "b"  || atom == "boron")         return 5; 
  if(atom == "c"  || atom == "carbon")        return 6; 
  if(atom == "n"  || atom == "nitrogen")      return 7; 
  if(atom == "o"  || atom == "oxygen")        return 8; 
  if(atom == "f"  || atom == "fluorine")      return 9; 
  if(atom == "ne" || atom == "neon")          return 10; 
  if(atom == "na" || atom == "sodium")        return 11; 
  if(atom == "mg" || atom == "magnesium")     return 12; 
  if(atom == "al" || atom == "aluminum")      return 13; 
  if(atom == "si" || atom == "silicon")       return 14; 
  if(atom == "p"  || atom == "phosphorus")    return 15; 
  if(atom == "s"  || atom == "sulfur")        return 16; 
  if(atom == "cl" || atom == "chlorine")      return 17; 
  if(atom == "ar" || atom == "argon")         return 18; 
  if(atom == "k"  || atom == "potassium")     return 19; 
  if(atom == "ca" || atom == "calcium")       return 20; 
  if(atom == "sc" || atom == "scandium")      return 21; 
  if(atom == "ti" || atom == "titanium")      return 22; 
  if(atom == "v"  || atom == "vanadium")      return 23; 
  if(atom == "cr" || atom == "chromium")      return 24; 
  if(atom == "mn" || atom == "manganese")     return 25; 
  if(atom == "fe" || atom == "iron")          return 26; 
  if(atom == "co" || atom == "cobalt")        return 27; 
  if(atom == "ni" || atom == "nickel")        return 28; 
  if(atom == "cu" || atom == "copper")        return 29; 
  if(atom == "zn" || atom == "zinc")          return 30; 
  if(atom == "ga" || atom == "gallium")       return 31; 
  if(atom == "ge" || atom == "germanium")     return 32; 
  if(atom == "as" || atom == "arsenic")       return 33; 
  if(atom == "se" || atom == "selenium")      return 34; 
  if(atom == "br" || atom == "bromine")       return 35; 
  if(atom == "kr" || atom == "krypton")       return 36; 
  if(atom == "rb" || atom == "rubidium")      return 37; 
  if(atom == "sr" || atom == "strontium")     return 38; 
  if(atom == "y"  || atom == "yttrium")       return 39; 
  if(atom == "zr" || atom == "zirconium")     return 40; 
  if(atom == "nb" || atom == "niobium")       return 41; 
  if(atom == "mo" || atom == "molybdenum")    return 42; 
  if(atom == "tc" || atom == "technetium")    return 43; 
  if(atom == "ru" || atom == "ruthenium")     return 44; 
  if(atom == "rh" || atom == "rhodium")       return 45; 
  if(atom == "pd" || atom == "palladium")     return 46; 
  if(atom == "ag" || atom == "silver")        return 47; 
  if(atom == "cd" || atom == "cadmium")       return 48; 
  if(atom == "in" || atom == "indium")        return 49; 
  if(atom == "sn" || atom == "tin")           return 50; 
  if(atom == "sb" || atom == "antimony")      return 51; 
  if(atom == "te" || atom == "tellurium")     return 52; 
  if(atom == "i"  || atom == "iodine")        return 53; 
  if(atom == "xe" || atom == "xenon")         return 54; 
  if(atom == "cs" || atom == "cesium")        return 55; 
  if(atom == "ba" || atom == "barium")        return 56; 
  if(atom == "la" || atom == "lanthanum")     return 57; 
  if(atom == "ce" || atom == "cerium")        return 58; 
  if(atom == "pr" || atom == "praseodymium")  return 59; 
  if(atom == "nd" || atom == "neodymium")     return 60; 
  if(atom == "pm" || atom == "promethium")    return 61; 
  if(atom == "sm" || atom == "samarium")      return 62; 
  if(atom == "eu" || atom == "europium")      return 63; 
  if(atom == "gd" || atom == "gadolinium")    return 64; 
  if(atom == "tb" || atom == "terbium")       return 65; 
  if(atom == "dy" || atom == "dysprosium")    return 66; 
  if(atom == "ho" || atom == "holmium")       return 67; 
  if(atom == "er" || atom == "erbium")        return 68; 
  if(atom == "tm" || atom == "thulium")       return 69; 
  if(atom == "yb" || atom == "ytterbium")     return 70; 
  if(atom == "lu" || atom == "lutetium")      return 71; 
  if(atom == "hf" || atom == "hafnium")       return 72; 
  if(atom == "ta" || atom == "tantalum")      return 73; 
  if(atom == "w"  || atom == "tungsten")      return 74; 
  if(atom == "re" || atom == "rhenium")       return 75; 
  if(atom == "os" || atom == "osmium")        return 76; 
  if(atom == "ir" || atom == "iridium")       return 77; 
  if(atom == "pt" || atom == "platinum")      return 78; 
  if(atom == "au" || atom == "gold")          return 79; 
  if(atom == "hg" || atom == "mercury")       return 80; 
  if(atom == "tl" || atom == "thallium")      return 81; 
  if(atom == "pb" || atom == "lead")          return 82; 
  if(atom == "bi" || atom == "bismuth")       return 83; 
  if(atom == "po" || atom == "polonium")      return 84; 
  if(atom == "at" || atom == "astatine")      return 85; 
  if(atom == "rn" || atom == "radon")         return 86; 
  if(atom == "fr" || atom == "francium")      return 87; 
  if(atom == "ra" || atom == "radium")        return 88; 
  if(atom == "ac" || atom == "actinium")      return 89; 
  if(atom == "th" || atom == "thorium")       return 90; 
  if(atom == "pa" || atom == "protactinium")  return 91; 
  if(atom == "u"  || atom == "uranium")       return 92; 
  if(atom == "np" || atom == "neptunium")     return 93; 
  if(atom == "pu" || atom == "plutonium")     return 94; 
  if(atom == "am" || atom == "americium")     return 95; 
  if(atom == "cm" || atom == "curium")        return 96; 
  if(atom == "bk" || atom == "berkelium")     return 97; 
  if(atom == "cf" || atom == "californium")   return 98; 
  if(atom == "es" || atom == "einsteinium")   return 99; 
  if(atom == "fm" || atom == "fermium")       return 100; 
  if(atom == "md" || atom == "mendelevium")   return 101; 
  if(atom == "no" || atom == "nobelium")      return 102; 
  if(atom == "lr" || atom == "lawrencium")    return 103; 
  if(atom == "rf" || atom == "rutherfordium") return 104; 
  if(atom == "db" || atom == "dubnium")       return 105; 
  if(atom == "sg" || atom == "seaborgium")    return 106; 
  if(atom == "bh" || atom == "bohrium")       return 107; 
  if(atom == "hs" || atom == "hassium")       return 108; 
  if(atom == "mt" || atom == "meitnerium")    return 109; 
  return 0;
}

void PrintOptions() {
  if(!verbose) return;
  
  cout << endl << endl;
  
  cout << "============================" << endl;
  cout << "Using the following options: " << endl;
  cout << "============================" << endl;
  
  if(invert) {
    cout << "WARNING: THE SELECTION WILL BE THE INVERSE OF:" << endl;
  }
  
  if(inputNeutrinoTypes.size() > 0) {
    cout << "Incoming neutrino types:" << endl;
    cout << "------------------------" << endl;
    for(vector<int>::const_iterator n = inputNeutrinoTypes.begin(); n != inputNeutrinoTypes.end(); ++n) {
      cout << *n << endl;
    }
    cout << "------------------------" << endl;
  }
  
  if(inputNeutrinoEnergyRanges.size() > 0) {
    cout << "Neutrino energy ranges:" << endl;
    cout << "-----------------------" << endl;
    for(vector<pair<double, double> >::const_iterator e = inputNeutrinoEnergyRanges.begin(); e != inputNeutrinoEnergyRanges.end(); ++e) {
      cout << "Min: " << e->first << " GeV, Max: " << e->second << " GeV" << endl;
    }
    cout << "-----------------------" << endl;
  }
  
  if(inputNuclearTargets.size() > 0) {
    cout << "Nuclear targets:" << endl;
    cout << "----------------" << endl;
    for(vector<int>::const_iterator t = inputNuclearTargets.begin(); t != inputNuclearTargets.end(); ++t) {
      cout << *t << endl;
    }
    cout << "----------------" << endl;
  }
  
  if(reactionCodes.size() > 0) {
    cout << "Neut-style reaction codes:" << endl;
    cout << "--------------------------" << endl;
    for(vector<int>::const_iterator r = reactionCodes.begin(); r != reactionCodes.end(); ++r) {
      cout << *r << endl;
    }
    cout << "--------------------------" << endl;
  }

  if(geoTargets.size() > 0) {
    cout << "Geometric target areas: " << endl;
    cout << "----------------------- " << endl;
    if(geoManager) {
      for(vector<string>::const_iterator g = geoTargets.begin(); g != geoTargets.end(); ++g) {
        cout << *g << endl;
      }
    } else {
      cout << "No geometry input!" << endl;
    }
    cout << "----------------------- " << endl;
    cout << endl;
  }

  if(outputParticles.size() > 0) {
    cout << "Output particle numbers [PDG:N]: " << endl;
    cout << "--------------------------------" << endl;
    if(countLeptons >= 0) cout << "Total number of outgoing leptons: " << countLeptons << endl;
    if(outputParticles.size() < 10) {
      for(map<int, int>::const_iterator i = outputParticles.begin(); i != outputParticles.end(); ++i) {
        cout << "[" << i->first << ":" << i->second << "], ";
      }
    }
    else {
      cout << outputParticles.size() << " output particles defined, will only enumarate non-zero ones:" << endl;
      for(map<int, int>::const_iterator i = outputParticles.begin(); i != outputParticles.end(); ++i) {
        if(i->second > 0) {
          cout << "[" << i->first << ":" << i->second << "], ";
        }
      }
    }
    cout << endl;
    cout << "--------------------------------" << endl;
  }
  
}