tutorials/block/Fun4All_G4_block.C

0001 #pragma once
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003 #include <fun4all/SubsysReco.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/Fun4AllInputManager.h>
0006 #include <fun4all/Fun4AllDstOutputManager.h>
0007 #include <fun4all/Fun4AllOutputManager.h>
0008 #include <fun4all/Fun4AllDummyInputManager.h>
0009 #include <g4detectors/PHG4BlockSubsystem.h>
0010 #include <g4eval/PHG4DSTReader.h>
0011 #include <g4main/PHG4ParticleGun.h>
0012 #include <g4main/PHG4Reco.h>
0013 #include <g4main/PHG4TruthSubsystem.h>
0014 #include <phool/recoConsts.h>
0015 R__LOAD_LIBRARY(libg4eval.so)
0016 R__LOAD_LIBRARY(libfun4all.so)
0017 R__LOAD_LIBRARY(libg4testbench.so)
0018 R__LOAD_LIBRARY(libg4detectors.so)
0019 #endif
0020
0021 int Fun4All_G4_block(const int nEvents = 10, const char *outfile=NULL)
0022 {
0023
0024   gSystem->Load("libfun4all");
0025   gSystem->Load("libg4detectors");
0026   gSystem->Load("libg4testbench");
0027   gSystem->Load("libg4eval");
0028   gSystem->Load("libg4histos");
0029
0030   ///////////////////////////////////////////
0031   // Make the Server
0032   //////////////////////////////////////////
0033   Fun4AllServer *se = Fun4AllServer::instance();
0034   //  se->Verbosity(1); // enable some blabbering
0035
0036   recoConsts *rc = recoConsts::instance();
0037 // uncomment and change number (or not)if you want to use a fixed seed
0038 //  rc->set_IntFlag("RANDOMSEED", 12345);
0039
0040   // particle gun
0041   PHG4ParticleGun *gun = new PHG4ParticleGun("PGUN");
0042   //  gun->set_name("anti_proton");
0043   gun->set_name("pi-");
0044   //  gun->set_name("proton");
0045   gun->set_vtx(0,0,0);
0046   gun->set_mom(0,0,1);
0047   se->registerSubsystem(gun);
0048
0049   // Fun4All G4 module
0050   PHG4Reco* g4Reco = new PHG4Reco();
0051   // no magnetic field
0052   g4Reco->set_field(0);
0053   // size of the world - every detector has to fit in here
0054   g4Reco->SetWorldSizeX(500);
0055   g4Reco->SetWorldSizeY(500);
0056   g4Reco->SetWorldSizeZ(2000);
0057   // shape of our world - it is a box
0058   g4Reco->SetWorldShape("G4BOX");
0059   // this is what our world is filled with
0060   g4Reco->SetWorldMaterial("G4_AIR");
0061   // Geant4 Physics list to use
0062   g4Reco->SetPhysicsList("QGSP_BERT");
0063
0064   // our block "detector", size is in cm
0065   double xsize = 200.;
0066   double ysize = 200.;
0067   double zsize = 400.;
0068   PHG4BlockSubsystem *box = new PHG4BlockSubsystem("box");
0069   box->set_double_param("size_x",xsize);
0070   box->set_double_param("size_y",ysize);
0071   box->set_double_param("size_z",zsize);
0072   box->set_double_param("place_z",zsize/2.+100);// shift box so we do not create particles in its center and shift by 10 so we can see the track of the incoming particle
0073   box->set_string_param("material","G4_POLYSTYRENE"); // material of box
0074   box->SetActive(); // it is an active volume - save G4Hits
0075   g4Reco->registerSubsystem(box);
0076
0077   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0078   g4Reco->registerSubsystem(truth);
0079
0080   se->registerSubsystem( g4Reco );
0081
0082   ///////////////////////////////////////////
0083   // Output
0084   ///////////////////////////////////////////
0085
0086   if (outfile)
0087     {
0088       Fun4AllOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",outfile);
0089       se->registerOutputManager(out);
0090
0091     }
0092   if (outfile)
0093     {
0094       // save a comprehensive  evaluation file
0095       PHG4DSTReader* ana = new PHG4DSTReader(
0096           string(outfile) + string("_DSTReader.root"));
0097       ana->set_save_particle(true);
0098       ana->set_load_all_particle(false);
0099       ana->set_load_active_particle(true);
0100       ana->set_save_vertex(true);
0101       if (nEvents > 0 && nEvents < 2)
0102         {
0103           ana->Verbosity(2);
0104         }
0105       ana->AddNode("box_0");
0106       se->registerSubsystem(ana);
0107     }
0108
0109    // input - we need a dummy to drive the event loop
0110   Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0111   se->registerInputManager( in );
0112
0113   // a quick evaluator to inspect on hit/particle/tower level
0114
0115
0116   if (nEvents > 0)
0117     {
0118       se->run(nEvents);
0119       // finish job - close and save output files
0120       se->End();
0121       std::cout << "All done" << std::endl;
0122
0123       // cleanup - delete the server and exit
0124       delete se;
0125       gSystem->Exit(0);
0126     }
0127    return 0;
0128
0129 }
0130
0131 PHG4ParticleGun *get_gun(const char *name="PGUN")
0132 {
0133   Fun4AllServer *se = Fun4AllServer::instance();
0134   PHG4ParticleGun *pgun = (PHG4ParticleGun  *) se->getSubsysReco(name);
0135   return pgun;
0136 }
0137