macros/tpc2019/Fun4All_G4_TPC.C

0001 #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 00, 0)
0002
0003 #include <fun4all/Fun4AllDstOutputManager.h>
0004 #include <fun4all/Fun4AllDummyInputManager.h>
0005 #include <fun4all/Fun4AllInputManager.h>
0006 #include <fun4all/Fun4AllOutputManager.h>
0007 #include <fun4all/Fun4AllServer.h>
0008 #include <fun4all/SubsysReco.h>
0009
0010 #include <fun4all/Fun4AllServer.h>
0011
0012 #include <g4eval/SvtxEvaluator.h>
0013 #include <g4eval/TrkrEvaluator.h>
0014
0015 #include <g4detectors/PHG4BlockSubsystem.h>
0016 #include <g4main/PHG4Reco.h>
0017 #include <phgeom/PHGeomFileImport.h>
0018
0019 #include <g4tpc/PHG4TpcDigitizer.h>
0020 #include <g4tpc/PHG4TpcElectronDrift.h>
0021 #include <g4tpc/PHG4TpcPadPlane.h>
0022 #include <g4tpc/PHG4TpcPadPlaneReadout.h>
0023 #include <g4tpc/PHG4TpcSubsystem.h>
0024 #include <trackreco/PHGenFitTrkFitter.h>
0025 #include <trackreco/PHGenFitTrkProp.h>
0026 #include <trackreco/PHHoughSeeding.h>
0027 #include <trackreco/PHInitVertexing.h>
0028 #include <trackreco/PHTrackSeeding.h>
0029 #include <trackreco/PHTruthTrackSeeding.h>
0030 #include <trackreco/PHTruthVertexing.h>
0031
0032 #include <g4eval/PHG4DSTReader.h>
0033
0034 #include <g4histos/G4HitNtuple.h>
0035
0036 #include <g4main/PHG4ParticleGenerator.h>
0037 #include <g4main/PHG4ParticleGun.h>
0038 #include <g4main/PHG4Reco.h>
0039 #include <g4main/PHG4SimpleEventGenerator.h>
0040 #include <g4main/PHG4TruthSubsystem.h>
0041
0042 #include <phool/recoConsts.h>
0043
0044 #include <tpc2019/TpcPrototypeClusterizer.h>
0045 #include <tpc2019/TpcPrototypeGenFitTrkFinder.h>
0046 #include <tpc2019/TpcPrototypeGenFitTrkFitter.h>
0047
0048 R__LOAD_LIBRARY(libcalo_reco.so)
0049 R__LOAD_LIBRARY(libfun4all.so)
0050 R__LOAD_LIBRARY(libg4tpc.so)
0051 R__LOAD_LIBRARY(libg4detectors.so)
0052 R__LOAD_LIBRARY(libg4eval.so)
0053 R__LOAD_LIBRARY(libg4histos.so)
0054 R__LOAD_LIBRARY(libg4testbench.so)
0055 R__LOAD_LIBRARY(libg4tpc.so)
0056 R__LOAD_LIBRARY(libg4intt.so)
0057 R__LOAD_LIBRARY(libg4mvtx.so)
0058 R__LOAD_LIBRARY(libg4hough.so)
0059 R__LOAD_LIBRARY(libg4eval.so)
0060 R__LOAD_LIBRARY(libintt.so)
0061 R__LOAD_LIBRARY(libmvtx.so)
0062 R__LOAD_LIBRARY(libtpc2019.so)
0063 R__LOAD_LIBRARY(libtrack_reco.so)
0064 #endif
0065
0066 int n_tpc_layer_inner = 0;
0067 int tpc_layer_rphi_count_inner = 1152;
0068 int n_tpc_layer_mid = 16;
0069 int n_tpc_layer_outer = 0;
0070 int n_gas_layer = n_tpc_layer_inner + n_tpc_layer_mid + n_tpc_layer_outer;
0071
0072 int Fun4All_G4_TPC(int nEvents = 100, bool eventDisp = false, int verbosity = 0)
0073 {
0074   gSystem->Load("libfun4all");
0075   gSystem->Load("libg4detectors");
0076   gSystem->Load("libg4testbench");
0077   gSystem->Load("libg4histos");
0078   gSystem->Load("libg4eval.so");
0079   gSystem->Load("libqa_modules");
0080   gSystem->Load("libg4tpc");
0081   gSystem->Load("libtrack_io.so");
0082   gSystem->Load("libfun4all.so");
0083   gSystem->Load("libg4detectors.so");
0084   gSystem->Load("libtpc2019.so");
0085   gSystem->Load("libg4eval.so");
0086   gSystem->Load("libfun4all.so");
0087   gSystem->Load("libg4detectors.so");
0088   gSystem->Load("libg4hough.so");
0089   gSystem->Load("libtrack_reco.so");
0090
0091   bool bh_on = false;  // the surrounding boxes need some further thinking
0092   bool dstreader = true;
0093   bool dstoutput = false;
0094
0095   const double TPCDriftLength = 40;
0096
0097   ///////////////////////////////////////////
0098   // Make the Server
0099   //////////////////////////////////////////
0100   Fun4AllServer *se = Fun4AllServer::instance();
0101   se->Verbosity(1);
0102   recoConsts *rc = recoConsts::instance();
0103   // only set this if you want a fixed random seed to make
0104   // results reproducible for testing
0105   rc->set_IntFlag("RANDOMSEED", 12345678);
0106
0107   // simulated setup sits at eta=1, theta=40.395 degrees
0108   double theta = 90 - 5;
0109   double phi = 180 + 360 / 12 / 2;
0110   // shift in x with respect to midrapidity setup
0111   double add_place_z = -TPCDriftLength * .5;
0112   // Test beam generator
0113   PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0114   gen->add_particles("proton", 1);  // mu-,e-,anti_proton,pi-
0115   gen->set_vertex_distribution_mean(0.0, 0.0, add_place_z);
0116   gen->set_vertex_distribution_width(0.0, .7, .7);  // Rough beam profile size @ 16 GeV measured by Abhisek
0117   gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
0118                                         PHG4SimpleEventGenerator::Gaus,
0119                                         PHG4SimpleEventGenerator::Gaus);  // Gauss beam profile
0120   double angle = theta * TMath::Pi() / 180.;
0121   double eta = -1. * TMath::Log(TMath::Tan(angle / 2.));
0122   gen->set_eta_range(eta - 0.001, eta + 0.001);                                          // 1mrad angular divergence
0123   gen->set_phi_range(TMath::Pi() * phi / 180 - 0.001, TMath::Pi() * phi / 180 + 0.001);  // 1mrad angular divergence
0124   const double momentum = 120;
0125   gen->set_p_range(momentum, momentum, momentum * 2e-2);  // 2% momentum smearing
0126   se->registerSubsystem(gen);
0127
0128   PHG4Reco *g4Reco = new PHG4Reco();
0129   g4Reco->set_field(0);
0130   //  g4Reco->SetPhysicsList("QGSP_BERT_HP"); // uncomment this line to enable the high-precision neutron simulation physics list, QGSP_BERT_HP
0131
0132   //----------------------------------------
0133   // TPC G4
0134   //----------------------------------------
0135
0136   PHG4TpcSubsystem *tpc = new PHG4TpcSubsystem("TPC");
0137   tpc->SetActive();
0138   tpc->SuperDetector("TPC");
0139   tpc->set_double_param("steplimits", 1);
0140   tpc->set_double_param("tpc_length", TPCDriftLength * 2);  // 2x 40 cm drift
0141   // By default uses "sPHENIX_TPC_Gas", defined in PHG4Reco. That is 90:10 Ne:C4
0142   tpc->SetAbsorberActive();
0143   g4Reco->registerSubsystem(tpc);
0144
0145   if (bh_on)
0146   {
0147     // BLACKHOLE, box surrounding the prototype to check for leakage
0148     PHG4BlockSubsystem *bh[5];
0149     // surrounding outer hcal
0150     // top
0151     bh[0] = new PHG4BlockSubsystem("bh1", 1);
0152     bh[0]->set_double_param("size_x", 270.);
0153     bh[0]->set_double_param("size_y", 0.01);
0154     bh[0]->set_double_param("size_z", 165.);
0155     bh[0]->set_double_param("place_x", 270. / 2.);
0156     bh[0]->set_double_param("place_y", 125. / 2.);
0157     // bottom
0158     bh[1] = new PHG4BlockSubsystem("bh2", 2);
0159     bh[1]->set_double_param("size_x", 270.);
0160     bh[1]->set_double_param("size_y", 0.01);
0161     bh[1]->set_double_param("size_z", 165.);
0162     bh[1]->set_double_param("place_x", 270. / 2.);
0163     bh[1]->set_double_param("place_y", -125. / 2.);
0164     // right side
0165     bh[2] = new PHG4BlockSubsystem("bh3", 3);
0166     bh[2]->set_double_param("size_x", 200.);
0167     bh[2]->set_double_param("size_y", 125.);
0168     bh[2]->set_double_param("size_z", 0.01);
0169     bh[2]->set_double_param("place_x", 200. / 2.);
0170     bh[2]->set_double_param("place_z", 165. / 2.);
0171     // left side
0172     bh[3] = new PHG4BlockSubsystem("bh4", 4);
0173     bh[3]->set_double_param("size_x", 270.);
0174     bh[3]->set_double_param("size_y", 125.);
0175     bh[3]->set_double_param("size_z", 0.01);
0176     bh[3]->set_double_param("place_x", 270. / 2.);
0177     bh[3]->set_double_param("place_z", -165. / 2.);
0178     // back
0179     bh[4] = new PHG4BlockSubsystem("bh5", 5);
0180     bh[4]->set_double_param("size_x", 0.01);
0181     bh[4]->set_double_param("size_y", 125.);
0182     bh[4]->set_double_param("size_z", 165.);
0183     bh[4]->set_double_param("place_x", 270.);
0184     for (int i = 0; i < 5; i++)
0185     {
0186       bh[i]->BlackHole();
0187       bh[i]->SetActive();
0188       bh[i]->SuperDetector("BlackHole");
0189       bh[i]->OverlapCheck(true);
0190       g4Reco->registerSubsystem(bh[i]);
0191     }
0192   }
0193   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0194   g4Reco->registerSubsystem(truth);
0195
0196   g4Reco->save_DST_geometry(false);
0197   se->registerSubsystem(g4Reco);
0198
0199   //swap out with the test beam geometry for the analysis stage
0200   PHGeomFileImport *import = new PHGeomFileImport("TpcPrototypeGeometry.gdml");
0201   se->registerSubsystem(import);
0202
0203   //=========================
0204   // setup Tpc readout for filling cells
0205   // g4tpc/PHG4TpcElectronDrift uses
0206   // g4tpc/PHG4TpcPadPlaneReadout
0207   //=========================
0208
0209   PHG4TpcPadPlane *padplane = new PHG4TpcPadPlaneReadout();
0210
0211   // The pad plane readout default is set in PHG4TpcPadPlaneReadout
0212   //only build mid-layer and to size
0213   padplane->set_int_param("tpc_minlayer_inner", 0);
0214   padplane->set_int_param("ntpc_layers_inner", 0);
0215   padplane->set_int_param("ntpc_layers_outer", 0);
0216   padplane->set_double_param("maxdriftlength", TPCDriftLength);
0217   padplane->set_int_param("ntpc_phibins_mid", 16 * 8 * 12);
0218
0219   PHG4TpcElectronDrift *edrift = new PHG4TpcElectronDrift();
0220   edrift->Detector("TPC");
0221   // fudge factors to get drphi 150 microns (in mid and outer Tpc) and dz 500 microns cluster resolution
0222   // They represent effects not due to ideal gas properties and ideal readout plane behavior
0223   // defaults are 0.12 and 0.15, they can be changed here to get a different resolution
0224   edrift->set_double_param("added_smear_trans", 0.12);
0225   edrift->set_double_param("added_smear_long", 0.15);
0226   edrift->registerPadPlane(padplane);
0227   edrift->Verbosity(verbosity);
0228   ;
0229   se->registerSubsystem(edrift);
0230
0231   // Tpc
0232   //====
0233   PHG4TpcDigitizer *digitpc = new PHG4TpcDigitizer();
0234   digitpc->SetTpcMinLayer(0);
0235   double ENC = 670.0;  // standard
0236   digitpc->SetENC(ENC);
0237   double ADC_threshold = 4.0 * ENC;
0238   digitpc->SetADCThreshold(ADC_threshold);  // 4 * ENC seems OK
0239
0240   cout << " Tpc digitizer: Setting ENC to " << ENC << " ADC threshold to " << ADC_threshold << endl;
0241
0242   se->registerSubsystem(digitpc);
0243
0244   // For the Tpc
0245   //==========
0246   TpcPrototypeClusterizer *tpcclusterizer = new TpcPrototypeClusterizer();
0247   tpcclusterizer->Verbosity(verbosity);
0248   ;
0249   se->registerSubsystem(tpcclusterizer);
0250
0251   //-------------
0252   // Tracking
0253   //------------
0254
0255   // This should be true for everything except testing wirh truth track seeding!
0256   const bool use_track_prop = true;
0257   if (use_track_prop)
0258   {
0259     // Find all clusters associated with each seed track
0260     TpcPrototypeGenFitTrkFinder *finder = new TpcPrototypeGenFitTrkFinder();
0261     finder->Verbosity(verbosity);
0262     finder->set_do_evt_display(eventDisp);
0263     finder->set_do_eval(true);
0264     se->registerSubsystem(finder);
0265   }
0266   else
0267   {
0268     //--------------------------------------------------
0269     // Track finding using truth information
0270     //--------------------------------------------------
0271
0272     PHInitVertexing *init_vtx = new PHTruthVertexing("PHTruthVertexing");
0273     init_vtx->Verbosity(verbosity);
0274     ;
0275     se->registerSubsystem(init_vtx);
0276
0277     // For each truth particle, create a track and associate clusters with it using truth information
0278     PHTruthTrackSeeding *pat_rec = new PHTruthTrackSeeding("PHTruthTrackSeeding");
0279     pat_rec->Verbosity(verbosity);
0280     ;
0281     pat_rec->set_min_clusters_per_track(10);
0282     se->registerSubsystem(pat_rec);
0283   }
0284   //
0285   //------------------------------------------------
0286   // Fitting of tracks using Kalman Filter
0287   //------------------------------------------------
0288
0289   TpcPrototypeGenFitTrkFitter *kalman = new TpcPrototypeGenFitTrkFitter();
0290   kalman->Verbosity(verbosity);
0291   kalman->set_do_evt_display(eventDisp);
0292   kalman->set_do_eval(true);
0293   se->registerSubsystem(kalman);
0294
0295   //----------------
0296   // Tracking evaluation
0297   //----------------
0298
0299   SvtxEvaluator *eval;
0300   eval = new SvtxEvaluator("SVTXEVALUATOR", "G4TPC_eval.root", "SvtxTrackMap", 0, 0, n_gas_layer);
0301   eval->do_cluster_eval(true);
0302   eval->do_g4hit_eval(true);
0303   eval->do_hit_eval(true);  // enable to see the hits that includes the chamber physics...
0304   eval->do_gpoint_eval(true);
0305   eval->do_eval_light(false);
0306   eval->scan_for_embedded(false);  // take all tracks if false - take only embedded tracks if true
0307   eval->Verbosity(verbosity);
0308   ;
0309   se->registerSubsystem(eval);
0310
0311   //----------------------
0312   // save a comprehensive  evaluation file
0313   //----------------------
0314   if (dstreader)
0315   {
0316     PHG4DSTReader *ana = new PHG4DSTReader(string("G4TPC_DSTReader.root"));
0317     ana->set_save_particle(true);
0318     ana->set_load_all_particle(true);
0319     ana->set_load_active_particle(true);
0320     ana->set_save_vertex(true);
0321
0322     ana->AddNode("ABSORBER_TPC");
0323     ana->AddNode("TPC");
0324     if (bh_on)
0325       ana->AddNode("BlackHole");  // add a G4Hit node
0326
0327     se->registerSubsystem(ana);
0328   }
0329
0330   if (dstoutput)
0331   {
0332     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", "G4TPC.root");
0333     se->registerOutputManager(out);
0334   }
0335
0336   Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
0337   se->registerInputManager(in);
0338   if (nEvents <= 0)
0339   {
0340     return 0;
0341   }
0342
0343   gSystem->ListLibraries();
0344   se->run(nEvents);
0345
0346   se->End();
0347
0348   //   std::cout << "All done" << std::endl;
0349   delete se;
0350   //   return 0;
0351   gSystem->Exit(0);
0352   return 0;
0353 }
0354
0355 // for using QuickTest to check if macro loads
0356 void RunLoadTest() {}