Fun4All_MCEventGen.C

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_MCEventGen.C

int Fun4All_MCEventGen(
                       const int nEvents = 100,
                       const char * outputFile = "G4MCEventGen.root"
                       )
{
  //===============
  // Input options
  //===============

  // read files in HepMC format (typically output from event generators like hijing or pythia)
  const bool readhepmc = false; // read HepMC files

  // read files in EICTree format generated by eicsmear package
  const bool readeictree = false;

  // Use particle generator Pythia 8
  const bool runpythia8 = false;

  // Use particle generator Pythia 6
  const bool runpythia6 = true;
  const char * pythia6configfile = "config_pythia6/phpythia6_ep.cfg";

  // Use particle generator HEPGen
  const bool runhepgen = false;

  // Use particle generator Sartre
  const bool runsartre = false;

  // Other options
  const bool do_dst_compress = true;

  // Option to save DST output file (for later use with Genat4 simulation)
  const bool do_DSTOutput = true;

  // Option to convert DST to human command readable TTree for quick poke around the outputs
  const bool do_DSTReader = false;

  // Option to save events in ASCII HepMC format
  const bool do_ASCIIOutput = false;


  //---------------
  // Load libraries
  //---------------
  gSystem->Load("libfun4all.so");
  gSystem->Load("libphhepmc.so");
  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4eval.so");

  //---------------
  // Fun4All server
  //---------------
  Fun4AllServer *se = Fun4AllServer::instance();
  se->Verbosity(0); // uncomment for batch production running with minimal output messages
  // se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running

  // just if we set some flags somewhere in this macro
  recoConsts *rc = recoConsts::instance();
  // By default every random number generator uses
  // PHRandomSeed() which reads /dev/urandom to get its seed
  // if the RANDOMSEED flag is set its value is taken as seed
  // You can either set this to a random value using PHRandomSeed()
  // which will make all seeds identical (not sure what the point of
  // this would be:
  //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
  // or set it to a fixed value so you can debug your code
  // rc->set_IntFlag("RANDOMSEED", 12345);

  /* Set world parameters in reco consts */
  rc->set_FloatFlag("WorldSizex", 1000.);
  rc->set_FloatFlag("WorldSizey", 1000.);
  rc->set_FloatFlag("WorldSizez", 1000.);
  rc->set_CharFlag("WorldShape", "G4Tubs");


  //-----------------
  // Event generation
  //-----------------

  if (readhepmc)
    {
    }
  else if (readeictree)
    {
      // this module is needed to read the EICTree style records into our G4 sims
      ReadEICFiles *eicr = new ReadEICFiles();
      eicr->OpenInputFile("data/eictree_milou_dvcs_10x250.root");

      se->registerSubsystem(eicr);
    }
  else if (runpythia8)
    {
      gSystem->Load("libPHPythia8.so");

      PHPythia8* pythia8 = new PHPythia8();
      // see coresoftware/generators/PHPythia8 for example config
      pythia8->set_config_file("phpythia8.cfg");
      se->registerSubsystem(pythia8);
    }
  else if (runpythia6)
    {
      gSystem->Load("libPHPythia6.so");

      PHPythia6 *pythia6 = new PHPythia6();
      pythia6->set_config_file( pythia6configfile );
      se->registerSubsystem(pythia6);
    }
  else if (runhepgen)
    {
      gSystem->Load("libsHEPGen.so");

      sHEPGen *hepgen = new sHEPGen();
      hepgen->set_datacard_file("config/hepgen_eic_dvcs.data");
      hepgen->set_momentum_electron(-10);
      hepgen->set_momentum_hadron(250);
      se->registerSubsystem(hepgen);
    }
  else if (runsartre)
    {
      // see coresoftware/generators/PHSartre/README for setup instructions
      // before running:
      // setenv SARTRE_DIR /opt/sphenix/core/sartre-1.20_root-5.34.36
      gSystem->Load("libPHSartre.so");

      PHSartre* mysartre = new PHSartre();
      // see coresoftware/generators/PHSartre for example config
      mysartre->set_config_file("config/sartre_ep.cfg");

      // particle trigger to enhance forward J/Psi -> ee
      //PHSartreParticleTrigger* pTrig = new PHSartreParticleTrigger("MySartreTrigger");
      //pTrig->AddParticles(-11);
      //pTrig->SetEtaHighLow(4.0,1.4);
      //pTrig->SetEtaHighLow(1.0,-1.1);  // central arm
      //pTrig->PrintConfig();
      //mysartre->register_trigger((PHSartreGenTrigger *)pTrig);
      se->registerSubsystem(mysartre);
    }

  /* Write DST output file */
  if ( do_DSTOutput )
    {
      Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
      se->registerOutputManager(out);
    }

  /* write DSTReader human readable output tree */
  if (do_DSTReader)
    {
      // load HepMCNodeReader
      HepMCNodeReader *hr = new HepMCNodeReader();
      se->registerSubsystem(hr);

      // load module to access truth particle information
      PHG4Reco* g4Reco = new PHG4Reco();
      PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
      g4Reco->registerSubsystem(truth);
      se->registerSubsystem( g4Reco );

      // save a comprehensive  evaluation file
      PHG4DSTReader* ana = new PHG4DSTReader(string(outputFile) + string("_DSTReader.root"));
      ana->set_save_particle(true);
      ana->set_load_all_particle(false);
      ana->set_load_active_particle(true);
      ana->set_save_vertex(true);

      se->registerSubsystem(ana);
    }

  /* Write HepMC ASCII output */
  else if ( do_ASCIIOutput )
    {
      Fun4AllHepMCOutputManager *asciiout = new Fun4AllHepMCOutputManager("HEPMCOUT",outputFile);
      se->registerOutputManager(asciiout);
    }

  //-----------------
  // Event processing
  //-----------------
  if (nEvents <= 0 && !readhepmc)
    {
      cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
      cout << "it will run forever, so I just return without running anything" << endl;
      return;
    }
  else
    {
      se->run(nEvents);

      se->End();
      std::cout << "All done" << std::endl;
      delete se;
      gSystem->Exit(0);
    }

}


void
G4Cmd(const char * cmd)
{
  Fun4AllServer *se = Fun4AllServer::instance();
  PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
  g4->ApplyCommand(cmd);
}