Fun4All_Generator_Display.C

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

#pragma once
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
#include <fun4all/SubsysReco.h>
#include <fun4all/Fun4AllServer.h>
#include <fun4all/Fun4AllInputManager.h>
#include <fun4all/Fun4AllDummyInputManager.h>
#include <g4main/PHG4ParticleGeneratorBase.h>
#include <g4main/PHG4ParticleGenerator.h>
#include <g4main/PHG4SimpleEventGenerator.h>
#include <g4main/PHG4ParticleGeneratorVectorMeson.h>
#include <g4main/PHG4ParticleGun.h>
#include <g4main/PHG4Reco.h>
#include <g4main/HepMCNodeReader.h>
#include <g4main/ReadEICFiles.h>
#include <phhepmc/Fun4AllHepMCInputManager.h>
#include <phool/recoConsts.h>
#include <phpythia6/PHPythia6.h>
#include <phpythia8/PHPythia8.h>
#include <phsartre/PHSartre.h>
R__LOAD_LIBRARY(libfun4all.so)
R__LOAD_LIBRARY(libg4testbench.so)
R__LOAD_LIBRARY(libPHPythia6.so)
R__LOAD_LIBRARY(libPHPythia8.so)
//R__LOAD_LIBRARY(libPHSartre.so)
#endif

using namespace std;

PHG4Reco *g4 = nullptr;

int Fun4All_Generator_Display(
  const int nEvents = 1,
  const std::string &inputFile = "input.root"
  )
{
  //===============
  // Input options
  //===============

  // Either:
  // read previously generated g4-hits files, in this case it opens a DST and skips
  // the simulations step completely. The G4Setup macro is only loaded to get information
  // about the number of layers used for the cell reco code
  //
  // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
  // Or:
  // read files in HepMC format (typically output from event generators like hijing or pythia)
  const bool readhepmc = false; // read HepMC files
  // Or:
  // read files in EICTree format generated by eicsmear package
  const bool readeictree = false;
  // Or:
  // Use Pythia 8
  const bool runpythia8 = false;
  // Or:
  // Use Pythia 6
  const bool runpythia6 = true;
  // Or:
  // Use Sartre - DO NOT USE RIGHT NOW
  const bool runsartre = false;




  // Besides the above flags. One can further choose to further put in following particles in Geant4 simulation
  // Use multi-particle generator (PHG4SimpleEventGenerator), see the code block below to choose particle species and kinematics
  const bool particles = false;
  // or gun/ very simple single particle gun generator
  const bool usegun = false;
  // Throw single Upsilons, may be embedded in Hijing by setting readhepmc flag also  (note, careful to set Z vertex equal to Hijing events)
  const bool upsilons = false;

  const double magfield = 1.5; // in T

  //---------------
  // Fun4All server
  //---------------

  Fun4AllServer *se = Fun4AllServer::instance();
//  se->Verbosity(1); // do some blabbering
  // 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
  // rc->set_IntFlag("RANDOMSEED", 12345);

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

  if (readeictree)
  {
    // this module is needed to read the EICTree style records into our G4 sims
    ReadEICFiles *eicr = new ReadEICFiles();
    eicr->OpenInputFile(inputFile);

    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();
    // see coresoftware/generators/PHPythia6 for example config
    pythia6->set_config_file("phpythia6.cfg");
    se->registerSubsystem(pythia6);
  }
  else if (runsartre)
  {
    // see coresoftware/generators/PHSartre/README for setup instructions
    gSystem->Load("libPHSartre.so");
    PHSartre* mysartre = new PHSartre();
    // see coresoftware/generators/PHSartre for example config
    mysartre->set_config_file("sartre.cfg");
    se->registerSubsystem(mysartre);
  }

  if(particles)
  {
    // toss low multiplicity dummy events
    PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
    gen->add_particles("pi-",1); // mu+,e+,proton,pi+,Upsilon
    //gen->add_particles("pi+",100); // 100 pion option
    if (readhepmc)
    {
      gen->set_reuse_existing_vertex(true);
      gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
    }
    else
    {
      gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
              PHG4SimpleEventGenerator::Uniform,
              PHG4SimpleEventGenerator::Uniform);
      gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
      gen->set_vertex_distribution_width(0.0, 0.0, 0.0);
    }
    gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
    gen->set_vertex_size_parameters(0.0, 0.0);
    gen->set_eta_range(-3, 3);
    gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
    gen->set_pt_range(0.1, 20.0);
    gen->Embed(1);
    gen->Verbosity(0);

    se->registerSubsystem(gen);
  }
  if (usegun)
  {
    // PHG4ParticleGun *gun = new PHG4ParticleGun();
    // gun->set_name("anti_proton");
    // gun->set_name("geantino");
    // gun->set_vtx(0, 0, 0);
    // gun->set_mom(10, 0, 0.01);
    // gun->AddParticle("geantino",1.7776,-0.4335,0.);
    // gun->AddParticle("geantino",1.7709,-0.4598,0.);
    // gun->AddParticle("geantino",2.5621,0.60964,0.);
    // gun->AddParticle("geantino",1.8121,0.253,0.);
    // se->registerSubsystem(gun);
    PHG4ParticleGenerator *pgen = new PHG4ParticleGenerator();
    pgen->set_name("mu-");
    pgen->set_z_range(0,0);
    pgen->set_eta_range(-4.0,0.0);
    pgen->set_mom_range(30,30);
    pgen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
    se->registerSubsystem(pgen);
  }

  // If "readhepMC" is also set, the Upsilons will be embedded in Hijing events, if 'particles" is set, the Upsilons will be embedded in whatever particles are thrown
  if(upsilons)
  {
    // run upsilons for momentum, dca performance, alone or embedded in Hijing

    PHG4ParticleGeneratorVectorMeson *vgen = new PHG4ParticleGeneratorVectorMeson();
    vgen->add_decay_particles("e+","e-",0); // i = decay id
    // event vertex
    if (readhepmc || particles)
    {
      vgen->set_reuse_existing_vertex(true);
    }

    // Note: this rapidity range completely fills the acceptance of eta = +/- 1 unit
    vgen->set_rapidity_range(-1.0, +1.0);
    vgen->set_pt_range(0.0, 10.0);

    int istate = 1;

    if(istate == 1)
    {
      // Upsilon(1S)
      vgen->set_mass(9.46);
      vgen->set_width(54.02e-6);
    }
    else if (istate == 2)
    {
      // Upsilon(2S)
      vgen->set_mass(10.0233);
      vgen->set_width(31.98e-6);
    }
    else
    {
      // Upsilon(3S)
      vgen->set_mass(10.3552);
      vgen->set_width(20.32e-6);
    }

    vgen->Verbosity(0);
    vgen->Embed(2);
    se->registerSubsystem(vgen);

    cout << "Upsilon generator for istate = " << istate << " created and registered "  << endl;
  }

  //--------------
  // Input Manager (HepMC or Dummy depending on input)
  //--------------

  if (readhepmc)
  {
    Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
    se->registerInputManager( in );
    se->fileopen( in->Name().c_str(), inputFile.c_str() );
  }
  else
  {
    // for single particle generators we just need something which drives
    // the event loop, the Dummy Input Mgr does just that
    Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
    se->registerInputManager( in );
  }
  // read-in HepMC events to Geant4 if there is any
  HepMCNodeReader *hr = new HepMCNodeReader();
  se->registerSubsystem(hr);
  g4 = new PHG4Reco();
  g4->set_rapidity_coverage(1.1); // according to drawings
  g4->set_field(magfield); 
  g4->save_DST_geometry(false); // saving the geometry crashes here
  se->registerSubsystem(g4);

// Start the display
  g4->InitRun(se->topNode());
  g4->ApplyDisplayAction();
  g4->ApplyCommand("/control/execute vis.mac");
// draw 1m long axis
  g4->ApplyCommand("/vis/scene/add/axes 0 0 0 100 cm");
  se->run(1);
// print some empty lines so the instructions stick out
  for (int i=0; i<5; i++)
  {
    cout << endl;
  }
  cout << "Type displaycmd() to see some commands to change the display" << endl;
  cout << "If you want to run more events, do:" << endl;
  cout << "Fun4AllServer *se = Fun4AllServer::instance();" << endl;
  cout << "se->run(1);" << endl;
  return 0;
}

void displaycmd()
{
  cout << "draw 1m axis: " << endl;
  cout << " g4->ApplyCommand(\"/vis/scene/add/axes 0 0 0 100 cm\")" << endl;
  cout << "zoom" << endl;
  cout << " g4->ApplyCommand(\"/vis/viewer/zoom 1\")" << endl;
  cout << "viewpoint:" << endl;
  cout << " g4->ApplyCommand(\"/vis/viewer/set/viewpointThetaPhi 0 0\")" << endl;
  cout << "panTo:" << endl;
  cout << " g4->ApplyCommand(\"/vis/viewer/panTo 0 0 cm\")" << endl;
  cout << "print to eps:" << endl;
  cout << " g4->ApplyCommand(\"/vis/ogl/printEPS\")" << endl;
  cout << "set background color:" << endl;
  cout << " g4->ApplyCommand(\"/vis/viewer/set/background white\")" << endl;
}