tutorials/IonGun/Fun4All_G4_IonGun.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 <g4histos/G4EdepNtuple.h>
0011 #include <g4histos/G4HitNtuple.h>
0012 #include <g4main/PHG4IonGun.h>
0013 #include <g4main/PHG4ParticleGun.h>
0014 #include <g4main/PHG4PrimaryGeneratorAction.h>
0015 #include <g4main/PHG4Reco.h>
0016 #include <g4main/PHG4TruthSubsystem.h>
0017 #include "GlobalVariables.h"
0018 #include "DisplayOn.C"
0019
0020 R__LOAD_LIBRARY(libfun4all.so)
0021 R__LOAD_LIBRARY(libg4testbench.so)
0022 R__LOAD_LIBRARY(libg4detectors.so)
0023 R__LOAD_LIBRARY(libg4histos)
0024 #endif
0025
0026 void Fun4All_G4_IonGun(int nEvents = 10)
0027 {
0028
0029   bool WriteDst = true; // set to true if yoy want to save everything in a Dst
0030
0031   ///////////////////////////////////////////
0032   // Make the Server
0033   //////////////////////////////////////////
0034   Fun4AllServer *se = Fun4AllServer::instance();
0035   //  se->Verbosity(1);
0036 // size of the box in cm
0037   double xsize = 20.;
0038   double ysize = 20.;
0039   double zsize = 10.;
0040 // needs to be registered before PHG4Reco, otherwise you generate
0041 // the ion after the G4 simulation is run leaving you with an empty detector
0042   igun = new PHG4IonGun();
0043   igun->SetA(197);
0044   igun->SetZ(79);
0045   igun->SetMom(0,0,1*197);
0046   igun->SetCharge(79);
0047   se->registerSubsystem(igun);
0048
0049   PHG4Reco* g4Reco = new PHG4Reco();
0050 // the default is a solenoidal field - we certainly do not want this
0051   g4Reco->set_field(0);
0052 // set the world size and shape, make it large enough for your volume to fit
0053 // but not much larger since G4 will track particles until they leave the
0054 // world volume (or they interact/decay)
0055   g4Reco->SetWorldSizeX(500);
0056   g4Reco->SetWorldSizeY(500);
0057   g4Reco->SetWorldSizeZ(500);
0058   g4Reco->SetWorldShape("G4BOX");
0059 // if you want vacuum use G4_Galactic (the thinnest material in G4)
0060   g4Reco->SetWorldMaterial("G4_AIR");
0061 // Choose an existing physics list. They are explicitely implemented in
0062 // our code, if you need another existing one, let us know.
0063 // BIC is binary intranuclear cascade, suitable for ions
0064   g4Reco->SetPhysicsList("QGSP_BIC");
0065
0066   PHG4BlockSubsystem *box = new PHG4BlockSubsystem("box1",1);
0067   box->set_double_param("size_x",xsize);
0068   box->set_double_param("size_y",ysize);
0069   box->set_double_param("size_z",zsize);
0070 // shift by zsize/2. puts ion gun just at the front of the target box
0071 // shifting by 10cm more gives some space to see the ion trail (but
0072 // you will incur energy loss in 10cm of air)
0073   box->set_double_param("place_z",zsize/2.+10);
0074   box->set_string_param("material","G4_WATER"); // material of target box
0075   // normally G4 determines the stepsize automatically according to what
0076 // physics process is chosen, leading to very coarse eloss resolution
0077 // this sets the step size to 1mm, caveat: if the step sizes are too small
0078 // the eloss calculation can be out of whack leading to very strange dEdx
0079   box->set_double_param("steplimits",0.1);
0080 // do not integrate energy loss in active vlomue, store each G4 step
0081 // as separate hit
0082   box->set_int_param("use_g4steps",1);
0083   box->SetActive();
0084   // this will check for overlaps during construction, if you have multiple
0085   // volumes in very close proximity you might want to run this once with
0086   // this flag enabled
0087   // box->OverlapCheck(1);
0088   // the name used to construct the hit node, used in the ntuple code
0089   // to identify the target volume (case sensitive)
0090   box->SuperDetector("box");
0091   g4Reco->registerSubsystem(box);
0092
0093   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0094   g4Reco->registerSubsystem(truth);
0095   se->registerSubsystem( g4Reco );
0096
0097
0098
0099 // first argument is the name, second the filename for the ntuple
0100   G4EdepNtuple *edn = new G4EdepNtuple("EdepNtuple","G4EdepNtuple.root");
0101   edn->AddNode("box",1); // connects hit node name with detid in ntuple (G4HIT_box is detid=1)
0102   se->registerSubsystem(edn);
0103   G4HitNtuple *hit = new G4HitNtuple("HitNtuple","G4HitNtuple.root");
0104   hit->AddNode("box",1); // connects hit node name with detid in ntuple (G4HIT_box is detid=1)
0105   se->registerSubsystem(hit);
0106   ///////////////////////////////////////////
0107   // IOManagers...
0108   ///////////////////////////////////////////
0109
0110   if (WriteDst)
0111   {
0112     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT","G4.root");
0113     se->registerOutputManager(out);
0114   }
0115   Fun4AllInputManager *in = new Fun4AllDummyInputManager( "DUMMY");
0116   se->registerInputManager( in );
0117 // only run if number of events is positive, otherwise quit here
0118 // so one can start the display and then run events by hand
0119   if (nEvents > 0)
0120   {
0121     se->run(nEvents);
0122     se->End();
0123     std::cout << "All done" << std::endl;
0124     delete se;
0125     gSystem->Exit(0);
0126   }
0127 }
0128
0129 PHG4IonGun *getgun()
0130 {
0131   return igun;
0132 }