tutorials/Momentum/Fun4All_G4_Momentum_Projection_Detectors.C

0001 #ifndef FUN4ALL_G4_MOMENTUM_PROJECTION_DETECTORS_C
0002 #define FUN4ALL_G4_MOMENTUM_PROJECTION_DETECTORS_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_Detectors(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, 1.5);
0045   gen->set_mom_range(2, 2);                          // GeV/c
0046   gen->set_phi_range(0., 90. / 180. * TMath::Pi());  // 0-90 deg
0047   se->registerSubsystem(gen);
0048
0049   PHG4Reco *g4Reco = new PHG4Reco();
0050   g4Reco->set_field(1.5);  // 1.5 T solenoidal field
0051
0052   double si_thickness[6] = {0.02, 0.02, 0.0625, 0.032, 0.032, 0.032};
0053   double svxrad[6] = {2.71, 4.63, 11.765, 25.46, 41.38, 63.66};
0054   double length[6] = {20., 20., 36., -1., -1., -1.};  // -1 use eta coverage to determine length
0055   PHG4CylinderSubsystem *cyl;
0056   // here is our silicon:
0057   for (int ilayer = 0; ilayer < 6; ilayer++)
0058   {
0059     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0060     cyl->set_double_param("radius", svxrad[ilayer]);
0061     cyl->set_string_param("material", "G4_Si");
0062     cyl->set_double_param("thickness", si_thickness[ilayer]);
0063     cyl->SetActive();
0064     cyl->SuperDetector("SVTX");
0065     if (length[ilayer] > 0)
0066     {
0067       cyl->set_double_param("length", length[ilayer]);
0068     }
0069     g4Reco->registerSubsystem(cyl);
0070   }
0071
0072   // projection volumes (thin cylinders)
0073   // if no material is given, world material is chosen (typically G4_AIR)
0074   double mycylinder1_radius = 2.;  // in cm
0075   cyl = new PHG4CylinderSubsystem("MyCylinder1");
0076   cyl->set_double_param("radius", mycylinder1_radius);
0077   cyl->set_double_param("thickness", 0.01);  // does not matter (but > 0)
0078   cyl->SuperDetector("MyCylinder1");
0079   cyl->set_double_param("length", 90.);
0080   cyl->SetActive();
0081   cyl->SaveAllHits();  // save all hits, also zero energy hits (which are normally discarded)
0082   g4Reco->registerSubsystem(cyl);
0083
0084   double mycylinder2_radius = 70.;  // in cm
0085   cyl = new PHG4CylinderSubsystem("MyCylinder2");
0086   cyl->set_double_param("radius", mycylinder2_radius);
0087   cyl->set_double_param("thickness", 0.01);  // does not matter (but > 0)
0088   cyl->SuperDetector("MyCylinder2");
0089   cyl->set_double_param("length", 90.);
0090   cyl->SetActive();
0091   cyl->SaveAllHits();
0092   g4Reco->registerSubsystem(cyl);
0093
0094   double myplane1_z = 100.;
0095   cyl = new PHG4CylinderSubsystem("MyPlane1");
0096   cyl->set_double_param("radius", 0);                    // 80 cm
0097   cyl->set_double_param("thickness", 75);                // for a cylindrical plane thickness is diameter
0098   cyl->set_double_param("length", 0.01);                 // for a cylindrical plane length is depth
0099   cyl->set_double_param("place_z", myplane1_z + 0.005);  // position in z, 1/2 depth needs to be added to put front at exactly 100 cm
0100   cyl->SuperDetector("MyPlane1");
0101   cyl->SetActive();
0102   cyl->SaveAllHits();
0103   g4Reco->registerSubsystem(cyl);
0104
0105   // Black hole swallows everything - prevent loopers from returning
0106   // to inner detectors
0107   cyl = new PHG4CylinderSubsystem("BlackHole", 0);
0108   cyl->set_double_param("radius", 80);      // 80 cm
0109   cyl->set_double_param("thickness", 0.1);  // does not matter (but > 0)
0110   cyl->SetActive();
0111   cyl->BlackHole();  // eats everything
0112   g4Reco->registerSubsystem(cyl);
0113
0114   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0115   g4Reco->registerSubsystem(truth);
0116
0117   se->registerSubsystem(g4Reco);
0118
0119   //---------------------------
0120   // fast pattern recognition and full Kalman filter
0121   // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
0122   //---------------------------
0123   PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
0124   kalman->set_use_vertex_in_fitting(false);
0125   kalman->set_sub_top_node_name("SVTX");
0126   kalman->set_trackmap_out_name("SvtxTrackMap");
0127
0128   //  add Si Trtacker
0129   kalman->add_phg4hits(
0130       "G4HIT_SVTX",                //      const std::string& phg4hitsNames,
0131       PHG4TrackFastSim::Cylinder,  //      const DETECTOR_TYPE phg4dettype,
0132       300e-4,                      //       radial-resolution [cm]
0133       30e-4,                       //        azimuthal-resolution [cm]
0134       1,                           //      z-resolution [cm]
0135       1,                           //      efficiency,
0136       0                            //      noise hits
0137   );
0138
0139   kalman->add_cylinder_state("MyCylinder1", mycylinder1_radius);
0140   kalman->add_cylinder_state("MyCylinder2", mycylinder2_radius);
0141   kalman->add_zplane_state("MyPlane1", myplane1_z);
0142
0143   se->registerSubsystem(kalman);
0144
0145   PHG4TrackFastSimEval *fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval");
0146   fast_sim_eval->set_filename(evalfile);
0147   fast_sim_eval->AddProjection("MyCylinder1");
0148   fast_sim_eval->AddProjection("MyCylinder2");
0149   fast_sim_eval->AddProjection("MyPlane1");
0150
0151   se->registerSubsystem(fast_sim_eval);
0152   //---------------------------
0153
0154   //---------------------------
0155   // output DST file for further offlien analysis
0156   //---------------------------
0157   if (!outfile.empty())
0158   {
0159     Fun4AllOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outfile);
0160     se->registerOutputManager(out);
0161   }
0162   Fun4AllInputManager *in = new Fun4AllDummyInputManager("ZEN");
0163   se->registerInputManager(in);
0164
0165   if (nEvents > 0)
0166   {
0167     se->run(nEvents);
0168     // finish job - close and save output files
0169     se->End();
0170     std::cout << "All done" << std::endl;
0171
0172     // cleanup - delete the server and exit
0173     delete se;
0174     gSystem->Exit(0);
0175   }
0176   return 0;
0177 }
0178
0179 PHG4ParticleGenerator *get_gen(const char *name = "PGENERATOR")
0180 {
0181   Fun4AllServer *se = Fun4AllServer::instance();
0182   PHG4ParticleGenerator *pgun = (PHG4ParticleGenerator *) se->getSubsysReco(name);
0183   return pgun;
0184 }
0185
0186 #endif