tutorials/Momentum/Fun4All_G4_Momentum_Projection.C

0001 #ifndef FUN4ALL_G4_MOMENTUM_PROJECTION_C
0002 #define FUN4ALL_G4_MOMENTUM_PROJECTION_C
0003
0004 #include <g4detectors/PHG4CylinderSubsystem.h>
0005
0006 #include <g4trackfastsim/PHG4TrackFastSim.h>
0007 #include <g4trackfastsim/PHG4TrackFastSimEval.h>
0008
0009 #include <g4main/PHG4ParticleGenerator.h>
0010 #include <g4main/PHG4ParticleGun.h>
0011 #include <g4main/PHG4Reco.h>
0012 #include <g4main/PHG4TruthSubsystem.h>
0013
0014 #include <fun4all/Fun4AllDstOutputManager.h>
0015 #include <fun4all/Fun4AllDummyInputManager.h>
0016 #include <fun4all/Fun4AllInputManager.h>
0017 #include <fun4all/Fun4AllOutputManager.h>
0018 #include <fun4all/Fun4AllServer.h>
0019 #include <fun4all/SubsysReco.h>
0020
0021 #include <phool/recoConsts.h>
0022
0023 R__LOAD_LIBRARY(libfun4all.so)
0024 R__LOAD_LIBRARY(libg4testbench.so)
0025 R__LOAD_LIBRARY(libg4detectors.so)
0026 R__LOAD_LIBRARY(libg4trackfastsim.so)
0027
0028 int Fun4All_G4_Momentum_Projection(const int nEvents = 1000, const string &evalfile = "FastTrackingEval.root", const string &outfile = "")
0029 {
0030   ///////////////////////////////////////////
0031   // Make the Server
0032   //////////////////////////////////////////
0033   Fun4AllServer *se = Fun4AllServer::instance();
0034   se->Verbosity(0);
0035
0036   recoConsts *rc = recoConsts::instance();
0037   // if you want to use a fixed seed for reproducible results
0038   //rc->set_IntFlag("RANDOMSEED", 12345); // if you want to use a fixed seed
0039   // PHG4ParticleGenerator generates particle
0040   // distributions in eta/phi/mom range
0041   PHG4ParticleGenerator *gen = new PHG4ParticleGenerator("PGENERATOR");
0042   gen->set_name("pi-");
0043   gen->set_vtx(0, 0, 0);
0044   gen->set_eta_range(0.5, 0.5);
0045   gen->set_mom_range(2, 2);                   // GeV/c
0046   gen->set_phi_range(0., 90. / 180. * M_PI);  // 0-90 deg
0047
0048   se->registerSubsystem(gen);
0049
0050   PHG4Reco *g4Reco = new PHG4Reco();
0051   g4Reco->set_field(1.5);  // 1.5 T solenoidal field
0052
0053   double si_thickness[6] = {0.02, 0.02, 0.0625, 0.032, 0.032, 0.032};
0054   double svxrad[6] = {2.71, 4.63, 11.765, 25.46, 41.38, 63.66};
0055   double length[6] = {20., 20., 36., -1., -1., -1.};  // -1 use eta coverage to determine length
0056   PHG4CylinderSubsystem *cyl;
0057   // here is our silicon:
0058   for (int ilayer = 0; ilayer < 6; ilayer++)
0059   {
0060     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0061     cyl->set_double_param("radius", svxrad[ilayer]);
0062     cyl->set_string_param("material", "G4_Si");
0063     cyl->set_double_param("thickness", si_thickness[ilayer]);
0064     cyl->SetActive();
0065     cyl->SuperDetector("SVTX");
0066     if (length[ilayer] > 0)
0067     {
0068       cyl->set_double_param("length", length[ilayer]);
0069     }
0070     g4Reco->registerSubsystem(cyl);
0071   }
0072
0073   // Black hole swallows everything - prevent loopers from returning
0074   // to inner detectors
0075   cyl = new PHG4CylinderSubsystem("BlackHole", 0);
0076   cyl->set_double_param("radius", 80);      // 80 cm
0077   cyl->set_double_param("thickness", 0.1);  // does not matter (but > 0)
0078   cyl->SetActive();
0079   cyl->BlackHole();  // eats everything
0080   g4Reco->registerSubsystem(cyl);
0081
0082   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0083   g4Reco->registerSubsystem(truth);
0084
0085   se->registerSubsystem(g4Reco);
0086
0087   //---------------------------
0088   // fast pattern recognition and full Kalman filter
0089   // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
0090   //---------------------------
0091   PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
0092   kalman->set_use_vertex_in_fitting(false);
0093   kalman->set_sub_top_node_name("SVTX");
0094   kalman->set_trackmap_out_name("SvtxTrackMap");
0095
0096   //  add Si Trtacker
0097   kalman->add_phg4hits(
0098       "G4HIT_SVTX",                //      const std::string& phg4hitsNames,
0099       PHG4TrackFastSim::Cylinder,  //      const DETECTOR_TYPE phg4dettype,
0100       300e-4,                      //       radial-resolution [cm]
0101       30e-4,                       //        azimuthal-resolution [cm]
0102       1,                           //      z-resolution [cm]
0103       1,                           //      efficiency,
0104       0                            //      noise hits
0105   );
0106
0107   kalman->add_cylinder_state("MyCylinder1", 2.);   // projection onto cylinder with radius = 2cm
0108   kalman->add_cylinder_state("MyCylinder2", 70.);  // projection onto cylinder with radius = 70cm
0109   kalman->add_zplane_state("MyPlane1", 100.);      // projection onto z-plane at 100cm
0110
0111   se->registerSubsystem(kalman);
0112
0113   PHG4TrackFastSimEval *fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval");
0114   fast_sim_eval->set_filename(evalfile);
0115   // Add the above projections to evaluation ntuple
0116   fast_sim_eval->AddProjection("MyCylinder1");
0117   fast_sim_eval->AddProjection("MyCylinder2");
0118   fast_sim_eval->AddProjection("MyPlane1");
0119
0120   se->registerSubsystem(fast_sim_eval);
0121   //---------------------------
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("JASMINE");
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