tutorials/Momentum/Fun4All_G4_Momentum_Projection_Calorimeters.C

0001 #ifndef FUN4ALL_G4_MOMENTUM_PROJECTION_CALORIMETERS_C
0002 #define FUN4ALL_G4_MOMENTUM_PROJECTION_CALORIMETERS_C
0003
0004 #include "G4_CEmc_Spacal.C"
0005
0006 #include <g4detectors/PHG4CylinderSubsystem.h>
0007
0008 #include <g4trackfastsim/PHG4TrackFastSim.h>
0009 #include <g4trackfastsim/PHG4TrackFastSimEval.h>
0010
0011 #include <g4main/PHG4ParticleGenerator.h>
0012 #include <g4main/PHG4ParticleGun.h>
0013 #include <g4main/PHG4Reco.h>
0014 #include <g4main/PHG4TruthSubsystem.h>
0015
0016 #include <fun4all/Fun4AllDstOutputManager.h>
0017 #include <fun4all/Fun4AllDummyInputManager.h>
0018 #include <fun4all/Fun4AllInputManager.h>
0019 #include <fun4all/Fun4AllOutputManager.h>
0020 #include <fun4all/Fun4AllServer.h>
0021 #include <fun4all/SubsysReco.h>
0022
0023 #include <phool/recoConsts.h>
0024
0025 R__LOAD_LIBRARY(libfun4all.so)
0026 R__LOAD_LIBRARY(libg4testbench.so)
0027 R__LOAD_LIBRARY(libg4detectors.so)
0028 R__LOAD_LIBRARY(libg4trackfastsim.so)
0029
0030 int Fun4All_G4_Momentum_Projection_Calorimeters(const int nEvents = 1000, const string &evalfile = "FastTrackingEval.root", const string &outfile = "")
0031 {
0032   ///////////////////////////////////////////
0033   // Make the Server
0034   //////////////////////////////////////////
0035   Fun4AllServer *se = Fun4AllServer::instance();
0036   se->Verbosity(0);
0037
0038   recoConsts *rc = recoConsts::instance();
0039   // if you want to use a fixed seed for reproducible results
0040   //  rc->set_IntFlag("RANDOMSEED", 12345); // if you want to use a fixed seed
0041   // PHG4ParticleGenerator generates particle
0042   // distributions in eta/phi/mom range
0043   PHG4ParticleGenerator *gen = new PHG4ParticleGenerator("PGENERATOR");
0044   gen->set_name("pi-");
0045   gen->set_vtx(0, 0, 0);
0046   gen->set_eta_range(-1, 1.);
0047   gen->set_mom_range(2, 2);                          // GeV/c
0048   gen->set_phi_range(0., 90. / 180. * TMath::Pi());  // 0-90 deg
0049   se->registerSubsystem(gen);
0050
0051   PHG4Reco *g4Reco = new PHG4Reco();
0052   g4Reco->set_field(1.5);  // 1.5 T solenoidal field
0053
0054   double si_thickness[6] = {0.02, 0.02, 0.0625, 0.032, 0.032, 0.032};
0055   double svxrad[6] = {2.71, 4.63, 11.765, 25.46, 41.38, 63.66};
0056   double length[6] = {20., 20., 36., -1., -1., -1.};  // -1 use eta coverage to determine length
0057   PHG4CylinderSubsystem *cyl;
0058   // here is our silicon:
0059   for (int ilayer = 0; ilayer < 6; ilayer++)
0060   {
0061     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0062     cyl->set_double_param("radius", svxrad[ilayer]);
0063     cyl->set_string_param("material", "G4_Si");
0064     cyl->set_double_param("thickness", si_thickness[ilayer]);
0065     cyl->SetActive();
0066     cyl->SuperDetector("SVTX");
0067     if (length[ilayer] > 0)
0068     {
0069       cyl->set_double_param("length", length[ilayer]);
0070     }
0071     g4Reco->registerSubsystem(cyl);
0072   }
0073
0074   CEmc(g4Reco, 25, 4);
0075
0076   // Black hole swallows everything - prevent loopers from returning
0077   // to inner detectors
0078   cyl = new PHG4CylinderSubsystem("BlackHole", 0);
0079   cyl->set_double_param("radius", 100);     // 80 cm
0080   cyl->set_double_param("thickness", 0.1);  // does not matter (but > 0)
0081   cyl->SetActive();
0082   cyl->BlackHole();  // eats everything
0083   g4Reco->registerSubsystem(cyl);
0084
0085   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0086   g4Reco->registerSubsystem(truth);
0087
0088   se->registerSubsystem(g4Reco);
0089
0090   CEMC_Cells();
0091   CEMC_Towers();
0092
0093   //---------------------------
0094   // fast pattern recognition and full Kalman filter
0095   // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
0096   //---------------------------
0097   PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
0098   kalman->set_use_vertex_in_fitting(false);
0099   kalman->set_sub_top_node_name("SVTX");
0100   kalman->set_trackmap_out_name("SvtxTrackMap");
0101
0102   //  add Si Trtacker
0103   kalman->add_phg4hits(
0104       "G4HIT_SVTX",                //      const std::string& phg4hitsNames,
0105       PHG4TrackFastSim::Cylinder,  //      const DETECTOR_TYPE phg4dettype,
0106       300e-4,                      //       radial-resolution [cm]
0107       30e-4,                       //        azimuthal-resolution [cm]
0108       1,                           //      z-resolution [cm]
0109       1,                           //      efficiency,
0110       0                            //      noise hits
0111   );
0112
0113   kalman->add_state_name("CEMC");
0114
0115   se->registerSubsystem(kalman);
0116
0117   PHG4TrackFastSimEval *fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval");
0118   fast_sim_eval->set_filename(evalfile);
0119   fast_sim_eval->AddProjection("CEMC");
0120
0121   se->registerSubsystem(fast_sim_eval);
0122
0123   //---------------------------
0124   // output DST file for further offlien analysis
0125   //---------------------------
0126   if (!outfile.empty())
0127   {
0128     Fun4AllOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outfile);
0129     se->registerOutputManager(out);
0130   }
0131   Fun4AllInputManager *in = new Fun4AllDummyInputManager("ANA");
0132   se->registerInputManager(in);
0133
0134   if (nEvents > 0)
0135   {
0136     se->run(nEvents);
0137     // finish job - close and save output files
0138     se->End();
0139     std::cout << "All done" << std::endl;
0140
0141     // cleanup - delete the server and exit
0142     delete se;
0143     gSystem->Exit(0);
0144   }
0145   return 0;
0146 }
0147
0148 PHG4ParticleGenerator *get_gen(const char *name = "PGENERATOR")
0149 {
0150   Fun4AllServer *se = Fun4AllServer::instance();
0151   PHG4ParticleGenerator *pgun = (PHG4ParticleGenerator *) se->getSubsysReco(name);
0152   return pgun;
0153 }
0154
0155 #endif