macros/tpc2019/ConstructGeometry.C

0001 // $Id: $
0002
0003 /*!
0004  * \file ConstructGeometry.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include <fun4all/Fun4AllDstOutputManager.h>
0012 #include <fun4all/Fun4AllDummyInputManager.h>
0013 #include <fun4all/Fun4AllInputManager.h>
0014 #include <fun4all/Fun4AllOutputManager.h>
0015 #include <fun4all/Fun4AllServer.h>
0016 #include <fun4all/SubsysReco.h>
0017
0018 #include <fun4all/Fun4AllServer.h>
0019
0020 #include <g4eval/SvtxEvaluator.h>
0021 #include <g4eval/TrkrEvaluator.h>
0022
0023 #include <g4detectors/PHG4BlockSubsystem.h>
0024 #include <g4detectors/PHG4CylinderSubsystem.h>
0025 #include <g4main/PHG4Reco.h>
0026
0027 #include <g4eval/PHG4DSTReader.h>
0028
0029 #include <g4histos/G4HitNtuple.h>
0030
0031 #include <g4main/PHG4ParticleGenerator.h>
0032 #include <g4main/PHG4ParticleGun.h>
0033 #include <g4main/PHG4Reco.h>
0034 #include <g4main/PHG4SimpleEventGenerator.h>
0035 #include <g4main/PHG4TruthSubsystem.h>
0036
0037 #include <phgeom/PHGeomUtility.h>
0038 #include <phool/recoConsts.h>
0039
0040 #include <TMath.h>
0041 #include <TSystem.h>
0042
0043 R__LOAD_LIBRARY(libcalo_reco.so)
0044 R__LOAD_LIBRARY(libfun4all.so)
0045 R__LOAD_LIBRARY(libg4tpc.so)
0046 R__LOAD_LIBRARY(libg4detectors.so)
0047 R__LOAD_LIBRARY(libg4eval.so)
0048 R__LOAD_LIBRARY(libg4histos.so)
0049 R__LOAD_LIBRARY(libg4testbench.so)
0050 R__LOAD_LIBRARY(libg4tpc.so)
0051 R__LOAD_LIBRARY(libg4intt.so)
0052 R__LOAD_LIBRARY(libg4mvtx.so)
0053 R__LOAD_LIBRARY(libg4hough.so)
0054 R__LOAD_LIBRARY(libg4eval.so)
0055 R__LOAD_LIBRARY(libintt.so)
0056 R__LOAD_LIBRARY(libmvtx.so)
0057 R__LOAD_LIBRARY(libtpc2019.so)
0058 R__LOAD_LIBRARY(libtrack_reco.so)
0059
0060 void ConstructGeometry()
0061 {
0062   const double placementR = 0.5 * (40 + 60);
0063   const double rotaitonZ = TMath::Pi() + TMath::Pi() * 2 / 12. / 2.;
0064   //  const double rotaitonZ = 0;
0065   const double driftLength = 40;
0066   const double cageRadius = 20;
0067
0068   gSystem->Load("libfun4all");
0069   gSystem->Load("libg4detectors");
0070   gSystem->Load("libg4testbench");
0071   gSystem->Load("libg4histos");
0072   gSystem->Load("libg4eval.so");
0073   gSystem->Load("libqa_modules");
0074   gSystem->Load("libg4tpc");
0075   gSystem->Load("libtrack_io.so");
0076   gSystem->Load("libfun4all.so");
0077   gSystem->Load("libg4detectors.so");
0078   gSystem->Load("libtpc2019.so");
0079   gSystem->Load("libg4eval.so");
0080   gSystem->Load("libfun4all.so");
0081   gSystem->Load("libg4detectors.so");
0082   gSystem->Load("libg4hough.so");
0083   gSystem->Load("libtrack_reco.so");
0084
0085   ///////////////////////////////////////////
0086   // Make the Server
0087   //////////////////////////////////////////
0088   Fun4AllServer *se = Fun4AllServer::instance();
0089   se->Verbosity(1);
0090   recoConsts *rc = recoConsts::instance();
0091   // only set this if you want a fixed random seed to make
0092   // results reproducible for testing
0093   // rc->set_IntFlag("RANDOMSEED",12345678);
0094
0095   // simulated setup sits at eta=1, theta=40.395 degrees
0096   double theta = 90;
0097   double phi = 180 + 360 / 12 / 2;
0098   // shift in x with respect to midrapidity setup
0099   double add_place_z = -driftLength * .5;
0100   // Test beam generator
0101   PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0102   gen->add_particles("proton", 1);  // mu-,e-,anti_proton,pi-
0103   gen->set_vertex_distribution_mean(0.0, 0.0, add_place_z);
0104   gen->set_vertex_distribution_width(0.0, .0, .0);  // Rough beam profile size @ 16 GeV measured by Abhisek
0105   gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
0106                                         PHG4SimpleEventGenerator::Gaus,
0107                                         PHG4SimpleEventGenerator::Gaus);  // Gauss beam profile
0108   double angle = theta * TMath::Pi() / 180.;
0109   double eta = -1. * TMath::Log(TMath::Tan(angle / 2.));
0110   gen->set_eta_range(eta - 0.001, eta + 0.001);                                          // 1mrad angular divergence
0111   gen->set_phi_range(TMath::Pi() * phi / 180 - 0.001, TMath::Pi() * phi / 180 + 0.001);  // 1mrad angular divergence
0112   const double momentum = 120;
0113   gen->set_p_range(momentum, momentum, momentum * 2e-2);  // 2% momentum smearing
0114   se->registerSubsystem(gen);
0115
0116   PHG4Reco *g4Reco = new PHG4Reco();
0117   g4Reco->set_field(0);
0118   //  g4Reco->SetPhysicsList("QGSP_BERT_HP"); // uncomment this line to enable the high-precision neutron simulation physics list, QGSP_BERT_HP
0119
0120   {
0121     PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("TPC_GasVol", 0);
0122     cyl->set_double_param("length", driftLength);
0123     cyl->set_double_param("place_x", placementR * TMath::Cos(rotaitonZ));
0124     cyl->set_double_param("place_y", placementR * TMath::Sin(rotaitonZ));
0125     cyl->set_double_param("place_z", -driftLength / 2);
0126     cyl->set_double_param("radius", 0.0);
0127     cyl->set_int_param("lengthviarapidity", 0);
0128     cyl->set_string_param("material", "sPHENIX_TPC_Gas");
0129     cyl->set_double_param("thickness", cageRadius);
0130     cyl->SuperDetector("TPC");
0131     cyl->SetActive();
0132     cyl->OverlapCheck(1);
0133     g4Reco->registerSubsystem(cyl);
0134   }
0135   {
0136     PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("TPC_FieldCage", 1);
0137     cyl->set_double_param("length", driftLength);
0138     cyl->set_double_param("place_x", placementR * TMath::Cos(rotaitonZ));
0139     cyl->set_double_param("place_y", placementR * TMath::Sin(rotaitonZ));
0140     cyl->set_double_param("place_z", -driftLength / 2);
0141     cyl->set_double_param("radius", cageRadius + 1e-4);
0142     cyl->set_int_param("lengthviarapidity", 0);
0143     cyl->set_string_param("material", "G4_Cu");
0144     cyl->set_double_param("thickness", 0.00347);
0145     cyl->SuperDetector("TPC_Support");
0146     cyl->SetActive();
0147     cyl->OverlapCheck(1);
0148     g4Reco->registerSubsystem(cyl);
0149   }
0150
0151   for (int sign = -1; sign <= 1; sign += 2)
0152   {
0153     const double endcap_thickness = 0.5;
0154
0155     PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("TPC_EndCap", 4 + sign);
0156     cyl->set_double_param("length", endcap_thickness);
0157     cyl->set_double_param("place_x", placementR * TMath::Cos(rotaitonZ));
0158     cyl->set_double_param("place_y", placementR * TMath::Sin(rotaitonZ));
0159     cyl->set_double_param("place_z", -driftLength / 2 + sign * ((driftLength / 2) + endcap_thickness / 2 + 1e-4));
0160     cyl->set_double_param("radius", 0.0);
0161     cyl->set_int_param("lengthviarapidity", 0);
0162     cyl->set_string_param("material", "G4_Al");
0163     cyl->set_double_param("thickness", cageRadius);
0164     cyl->SuperDetector("TPC_Support");
0165     cyl->SetActive();
0166     cyl->OverlapCheck(1);
0167     g4Reco->registerSubsystem(cyl);
0168   }
0169
0170   for (int i = 0; i < 3; ++i)
0171   {
0172     const double GEMLocation = 0;
0173     const double GEMSpacing = -20;
0174
0175     PHG4BlockSubsystem *cyl = new PHG4BlockSubsystem("GEM", i);
0176     cyl->set_double_param("size_x", 0.5);
0177     cyl->set_double_param("size_y", 40);
0178     cyl->set_double_param("size_z", 40);
0179     cyl->set_double_param("place_x", (GEMLocation + GEMSpacing * i) * TMath::Cos(rotaitonZ));
0180     cyl->set_double_param("place_y", (GEMLocation + GEMSpacing * i) * TMath::Sin(rotaitonZ));
0181     cyl->set_double_param("place_z", -driftLength / 2);
0182     cyl->set_double_param("rot_z", -rotaitonZ * 180 / TMath::Pi());
0183     cyl->set_int_param("lengthviarapidity", 0);
0184     cyl->set_string_param("material", "ePHEINX_TPC_Gas"); // ePHENIX TPC is Ar + CO2 gas
0185     cyl->SuperDetector("GEM");
0186     cyl->SetActive();
0187     cyl->OverlapCheck(1);
0188     g4Reco->registerSubsystem(cyl);
0189   }
0190
0191   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0192   g4Reco->registerSubsystem(truth);
0193
0194   se->registerSubsystem(g4Reco);
0195
0196   Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
0197   se->registerInputManager(in);
0198
0199   char cmd[100];
0200   g4Reco->InitRun(se->topNode());
0201   g4Reco->ApplyDisplayAction();
0202   sprintf(cmd, "/control/execute vis.mac");
0203   g4Reco->ApplyCommand(cmd);
0204
0205   se->run(1);
0206   g4Reco->ApplyCommand("/vis/scene/add/axes 0 0 0 50 cm");
0207   g4Reco->ApplyCommand("/vis/viewer/zoom 1");
0208
0209   PHGeomUtility::ExportGeomtry(se->topNode(), "TpcPrototypeGeometry.root");
0210   PHGeomUtility::ExportGeomtry(se->topNode(), "TpcPrototypeGeometry.gdml");
0211 }