tutorials/eventgenerator_display/Fun4All_Generator_Display.C

0001 #ifndef MACRO_FUN4ALLGENERATORDISPLAY_C
0002 #define MACRO_FUN4ALLGENERATORDISPLAY_C
0003
0004 #include <fun4all/SubsysReco.h>
0005 #include <fun4all/Fun4AllServer.h>
0006 #include <fun4all/Fun4AllInputManager.h>
0007 #include <fun4all/Fun4AllDummyInputManager.h>
0008 #include <g4main/PHG4ParticleGeneratorBase.h>
0009 #include <g4main/PHG4ParticleGenerator.h>
0010 #include <g4main/PHG4SimpleEventGenerator.h>
0011 #include <g4main/PHG4ParticleGeneratorVectorMeson.h>
0012 #include <g4main/PHG4ParticleGun.h>
0013 #include <g4main/CosmicSpray.h>
0014 #include <g4main/PHG4Reco.h>
0015 #include <g4main/HepMCNodeReader.h>
0016 #include <g4main/ReadEICFiles.h>
0017 #include <phhepmc/Fun4AllHepMCInputManager.h>
0018 #include <phool/recoConsts.h>
0019 #include <phpythia6/PHPythia6.h>
0020 #include <phpythia8/PHPythia8.h>
0021 #include <phsartre/PHSartre.h>
0022
0023 R__LOAD_LIBRARY(libfun4all.so)
0024 R__LOAD_LIBRARY(libg4testbench.so)
0025 R__LOAD_LIBRARY(libPHPythia6.so)
0026 R__LOAD_LIBRARY(libPHPythia8.so)
0027 R__LOAD_LIBRARY(libPHSartre.so)
0028
0029 PHG4Reco *g4 = nullptr;
0030
0031 int Fun4All_Generator_Display(
0032   const int nEvents = 1,
0033   const std::string &inputFile = "input.root"
0034   )
0035 {
0036   //===============
0037   // Input options
0038   //===============
0039
0040   // Either:
0041   // read previously generated g4-hits files, in this case it opens a DST and skips
0042   // the simulations step completely. The G4Setup macro is only loaded to get information
0043   // about the number of layers used for the cell reco code
0044   //
0045   // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
0046   // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0047   // Or:
0048   // read files in HepMC format (typically output from event generators like hijing or pythia)
0049   const bool readhepmc = false; // read HepMC files
0050   // Or:
0051   // read files in EICTree format generated by eicsmear package
0052   const bool readeictree = false;
0053   // Or:
0054   // Use Pythia 8
0055   const bool runpythia8 = true;
0056   // Or:
0057   // Use Pythia 6
0058   const bool runpythia6 = false;
0059   // Or:
0060   // Use Sartre - DO NOT USE RIGHT NOW
0061   const bool runsartre = false;
0062
0063
0064
0065
0066   // Besides the above flags. One can further choose to further put in following particles in Geant4 simulation
0067   // Use multi-particle generator (PHG4SimpleEventGenerator), see the code block below to choose particle species and kinematics
0068   const bool particles = false;
0069   // or gun/ very simple single particle gun generator
0070   const bool usegun = false;
0071   // Throw single Upsilons, may be embedded in Hijing by setting readhepmc flag also  (note, careful to set Z vertex equal to Hijing events)
0072   const bool upsilons = false;
0073
0074   const bool cosmic = false;
0075
0076   const double magfield = 1.4; // in T
0077
0078   //---------------
0079   // Fun4All server
0080   //---------------
0081
0082   Fun4AllServer *se = Fun4AllServer::instance();
0083 //  se->Verbosity(1); // do some blabbering
0084   // just if we set some flags somewhere in this macro
0085   recoConsts *rc = recoConsts::instance();
0086   // By default every random number generator uses
0087   // PHRandomSeed() which reads /dev/urandom to get its seed
0088   // if the RANDOMSEED flag is set its value is taken as seed
0089   // rc->set_IntFlag("RANDOMSEED", 12345);
0090
0091   //-----------------
0092   // Event generation
0093   //-----------------
0094
0095   if (readeictree)
0096   {
0097     // this module is needed to read the EICTree style records into our G4 sims
0098     ReadEICFiles *eicr = new ReadEICFiles();
0099     eicr->OpenInputFile(inputFile);
0100
0101     se->registerSubsystem(eicr);
0102   }
0103   else if (runpythia8)
0104   {
0105     PHPythia8* pythia8 = new PHPythia8();
0106     // see coresoftware/generators/PHPythia8 for example config
0107     pythia8->set_config_file("phpythia8.cfg");
0108     se->registerSubsystem(pythia8);
0109   }
0110   else if (runpythia6)
0111   {
0112     PHPythia6 *pythia6 = new PHPythia6();
0113     // see coresoftware/generators/PHPythia6 for example config
0114     pythia6->set_config_file("phpythia6.cfg");
0115     se->registerSubsystem(pythia6);
0116   }
0117   else if (runsartre)
0118   {
0119     // see coresoftware/generators/PHSartre/README for setup instructions
0120     PHSartre* mysartre = new PHSartre();
0121     // see coresoftware/generators/PHSartre for example config
0122     mysartre->set_config_file("sartre.cfg");
0123     se->registerSubsystem(mysartre);
0124   }
0125
0126   if(particles)
0127   {
0128     // toss low multiplicity dummy events
0129     PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0130     gen->add_particles("pi-",1); // mu+,e+,proton,pi+,Upsilon
0131     //gen->add_particles("pi+",100); // 100 pion option
0132     if (readhepmc)
0133     {
0134       gen->set_reuse_existing_vertex(true);
0135       gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0136     }
0137     else
0138     {
0139       gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0140                         PHG4SimpleEventGenerator::Uniform,
0141                         PHG4SimpleEventGenerator::Uniform);
0142       gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
0143       gen->set_vertex_distribution_width(0.0, 0.0, 0.0);
0144     }
0145     gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0146     gen->set_vertex_size_parameters(0.0, 0.0);
0147     gen->set_eta_range(-3, 3);
0148     gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0149     gen->set_pt_range(0.1, 20.0);
0150     gen->Embed(1);
0151     gen->Verbosity(0);
0152
0153     se->registerSubsystem(gen);
0154   }
0155   if (usegun)
0156   {
0157     // PHG4ParticleGun *gun = new PHG4ParticleGun();
0158     // gun->set_name("anti_proton");
0159     // gun->set_name("geantino");
0160     // gun->set_vtx(0, 0, 0);
0161     // gun->set_mom(10, 0, 0.01);
0162     // gun->AddParticle("geantino",1.7776,-0.4335,0.);
0163     // gun->AddParticle("geantino",1.7709,-0.4598,0.);
0164     // gun->AddParticle("geantino",2.5621,0.60964,0.);
0165     // gun->AddParticle("geantino",1.8121,0.253,0.);
0166     // se->registerSubsystem(gun);
0167     PHG4ParticleGenerator *pgen = new PHG4ParticleGenerator();
0168     pgen->set_name("mu-");
0169     pgen->set_z_range(0,0);
0170     pgen->set_eta_range(-4.0,0.0);
0171     pgen->set_mom_range(30,30);
0172     pgen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0173     se->registerSubsystem(pgen);
0174   }
0175   if (cosmic)
0176   {
0177     CosmicSpray *cosmo = new CosmicSpray("CosmicGun",500.);
0178     se->registerSubsystem(cosmo);
0179   }
0180
0181   // 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
0182   if(upsilons)
0183   {
0184     // run upsilons for momentum, dca performance, alone or embedded in Hijing
0185
0186     PHG4ParticleGeneratorVectorMeson *vgen = new PHG4ParticleGeneratorVectorMeson();
0187     vgen->add_decay_particles("e+","e-",0); // i = decay id
0188     // event vertex
0189     if (readhepmc || particles)
0190     {
0191       vgen->set_reuse_existing_vertex(true);
0192     }
0193
0194     // Note: this rapidity range completely fills the acceptance of eta = +/- 1 unit
0195     vgen->set_rapidity_range(-1.0, +1.0);
0196     vgen->set_pt_range(0.0, 10.0);
0197
0198     int istate = 1;
0199
0200     if(istate == 1)
0201     {
0202       // Upsilon(1S)
0203       vgen->set_mass(9.46);
0204       vgen->set_width(54.02e-6);
0205     }
0206     else if (istate == 2)
0207     {
0208       // Upsilon(2S)
0209       vgen->set_mass(10.0233);
0210       vgen->set_width(31.98e-6);
0211     }
0212     else
0213     {
0214       // Upsilon(3S)
0215       vgen->set_mass(10.3552);
0216       vgen->set_width(20.32e-6);
0217     }
0218
0219     vgen->Verbosity(0);
0220     vgen->Embed(2);
0221     se->registerSubsystem(vgen);
0222
0223     cout << "Upsilon generator for istate = " << istate << " created and registered "  << endl;
0224   }
0225
0226   //--------------
0227   // Input Manager (HepMC or Dummy depending on input)
0228   //--------------
0229
0230   if (readhepmc)
0231   {
0232     Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0233     se->registerInputManager( in );
0234     se->fileopen( in->Name(), inputFile );
0235   }
0236   else
0237   {
0238     // for single particle generators we just need something which drives
0239     // the event loop, the Dummy Input Mgr does just that
0240     Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0241     se->registerInputManager( in );
0242   }
0243   // read-in HepMC events to Geant4 if there is any
0244   HepMCNodeReader *hr = new HepMCNodeReader();
0245   se->registerSubsystem(hr);
0246
0247   g4 = new PHG4Reco();
0248   g4->set_rapidity_coverage(1.1); // according to drawings
0249   g4->set_field(magfield);
0250   g4->save_DST_geometry(false); // saving the geometry crashes here
0251   se->registerSubsystem(g4);
0252
0253 // Start the display
0254
0255   g4->InitRun(se->topNode());
0256   g4->ApplyDisplayAction();
0257   g4->ApplyCommand("/control/execute vis.mac");
0258 // draw 1m long axis
0259   g4->ApplyCommand("/vis/scene/add/axes 0 0 0 100 cm");
0260
0261   se->run(1);
0262 // print some empty lines so the instructions stick out
0263   for (int i=0; i<5; i++)
0264   {
0265     cout << endl;
0266   }
0267   cout << "Type displaycmd() to see some commands to change the display" << endl;
0268   cout << "If you want to run more events, do:" << endl;
0269   cout << "Fun4AllServer *se = Fun4AllServer::instance();" << endl;
0270   cout << "se->run(1);" << endl;
0271   return 0;
0272 }
0273
0274 void displaycmd()
0275 {
0276   cout << "draw 1m axis: " << endl;
0277   cout << " g4->ApplyCommand(\"/vis/scene/add/axes 0 0 0 100 cm\")" << endl;
0278   cout << "zoom" << endl;
0279   cout << " g4->ApplyCommand(\"/vis/viewer/zoom 1\")" << endl;
0280   cout << "viewpoint:" << endl;
0281   cout << " g4->ApplyCommand(\"/vis/viewer/set/viewpointThetaPhi 0 0\")" << endl;
0282   cout << "panTo:" << endl;
0283   cout << " g4->ApplyCommand(\"/vis/viewer/panTo 0 0 cm\")" << endl;
0284   cout << "print to eps:" << endl;
0285   cout << " g4->ApplyCommand(\"/vis/ogl/printEPS\")" << endl;
0286   cout << "set background color:" << endl;
0287   cout << " g4->ApplyCommand(\"/vis/viewer/set/background white\")" << endl;
0288 }
0289
0290 #endif