MvtxTelescope/macros/Fun4All_G4_LeadBrick.C

0001
0002 #include <iostream>
0003
0004 using namespace std;
0005
0006 int Fun4All_G4_LeadBrick(const int nEvents = 10, const char *outfile = NULL)
0007 {
0008   gSystem->Load("libfun4all");
0009   gSystem->Load("libg4detectors");
0010   gSystem->Load("libg4testbench");
0011   gSystem->Load("libg4eval");
0012
0013   ///////////////////////////////////////////
0014   // Make the Server
0015   //////////////////////////////////////////
0016   Fun4AllServer *se = Fun4AllServer::instance();
0017   se->Verbosity(1);
0018
0019   // particle gun
0020   PHG4ParticleGun *gun = new PHG4ParticleGun("PGUN");
0021   //  gun->set_name("anti_proton");
0022   gun->set_name("proton");
0023   //  gun->set_name("mu-");
0024   //  gun->set_name("proton");
0025   gun->set_vtx(0, 0, -100);
0026   gun->set_mom(0, 0, 120);
0027   se->registerSubsystem(gun);
0028
0029   // Fun4All G4 module
0030   PHG4Reco *g4Reco = new PHG4Reco();
0031   // no magnetic field
0032   g4Reco->set_field(0);
0033   // size of the world - every detector has to fit in here
0034   g4Reco->SetWorldSizeX(200);
0035   g4Reco->SetWorldSizeY(200);
0036   g4Reco->SetWorldSizeZ(200);
0037   // shape of our world - it is a box
0038   g4Reco->SetWorldShape("G4BOX");
0039   // this is what our world is filled with
0040   g4Reco->SetWorldMaterial("G4_AIR");
0041   // Geant4 Physics list to use
0042   g4Reco->SetPhysicsList("FTFP_BERT");
0043
0044   // our block "detector", size is in cm
0045   double xsize = 200.;
0046   double ysize = 200.;
0047   double zsize = 400.;
0048   PHG4BlockSubsystem *box = new PHG4BlockSubsystem("LeadBrick");
0049   box->set_double_param("size_x", 20);
0050   box->set_double_param("size_y", 20);
0051   box->set_double_param("size_z", 10);
0052   box->set_double_param("place_z", 0);         // 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
0053   box->set_string_param("material", "G4_Pb");  // material of box
0054   box->SetActive(0);                           // if it is an active volume - save G4Hits
0055   g4Reco->registerSubsystem(box);
0056
0057   for (int stave = 0; stave < 4; ++stave)
0058   {
0059     box = new PHG4BlockSubsystem("MVTX", stave);
0060     box->SuperDetector("MVTX");
0061     box->set_double_param("size_x", 2);
0062     box->set_double_param("size_y", 1);
0063     box->set_double_param("size_z", 50e-4);              // 50us
0064     box->set_double_param("place_z", 10 * (stave + 1));  // 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
0065     box->set_string_param("material", "G4_Si");          // material of box
0066     box->SetActive(1);                                   // it is an active volume - save G4Hits
0067     g4Reco->registerSubsystem(box);
0068   }
0069
0070   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0071   g4Reco->registerSubsystem(truth);
0072
0073   se->registerSubsystem(g4Reco);
0074
0075   ///////////////////////////////////////////
0076   // Output
0077   ///////////////////////////////////////////
0078
0079   // save a comprehensive  evaluation file
0080   PHG4DSTReader *ana = new PHG4DSTReader(
0081       string("LeadBrick_DSTReader.root"));
0082   ana->set_save_particle(true);
0083   ana->set_load_all_particle(false);
0084   ana->set_load_active_particle(true);
0085   ana->set_save_vertex(true);
0086   ana->AddNode("LeadBrick");
0087   ana->AddNode("MVTX");
0088   se->registerSubsystem(ana);
0089
0090   // input - we need a dummy to drive the event loop
0091   Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
0092   se->registerInputManager(in);
0093
0094   // a quick evaluator to inspect on hit/particle/tower level
0095
0096   if (nEvents > 0)
0097   {
0098     se->run(nEvents);
0099     // finish job - close and save output files
0100     se->End();
0101     std::cout << "All done" << std::endl;
0102
0103     // cleanup - delete the server and exit
0104     delete se;
0105     gSystem->Exit(0);
0106   }
0107   return;
0108 }
0109
0110 PHG4ParticleGun *get_gun(const char *name = "PGUN")
0111 {
0112   Fun4AllServer *se = Fun4AllServer::instance();
0113   PHG4ParticleGun *pgun = (PHG4ParticleGun *) se->getSubsysReco(name);
0114   return pgun;
0115 }