FastTrackingEval/macros/Fun4All_G4_fsPHENIX_FastSim.C

0001
0002 int Fun4All_G4_fsPHENIX_FastSim(
0003                const int nEvents = 2,
0004                const char * inputFile = "/gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.5/single_particle/spacal1d/fieldmap/G4Hits_sPHENIX_e-_eta0_16GeV.root",
0005                const char * outputFile = "DST_PHG4TrackFastSim.root"
0006                )
0007 {
0008   //===============
0009   // Input options
0010   //===============
0011
0012   // Either:
0013   // read previously generated g4-hits files, in this case it opens a DST and skips
0014   // the simulations step completely. The G4Setup macro is only loaded to get information
0015   // about the number of layers used for the cell reco code
0016   const bool readhits = false;
0017   // Or:
0018   // read files in HepMC format (typically output from event generators like hijing or pythia)
0019   const bool readhepmc = false; // read HepMC files
0020   // Or:
0021   // Use particle generator
0022   const bool runpythia8 = false;
0023   const bool runpythia6 = false;
0024
0025   //======================
0026   // What to run
0027   //======================
0028
0029   bool do_bbc = false;
0030
0031   bool do_pipe = false;
0032
0033   bool do_svtx = false;
0034   bool do_svtx_cell = false;
0035   bool do_svtx_track = false;
0036   bool do_svtx_eval = false;
0037
0038   bool do_preshower = false;
0039
0040   bool do_cemc = false;
0041   bool do_cemc_cell = false;
0042   bool do_cemc_twr = false;
0043   bool do_cemc_cluster = false;
0044   bool do_cemc_eval = false;
0045
0046   bool do_hcalin = false;
0047   bool do_hcalin_cell = false;
0048   bool do_hcalin_twr = false;
0049   bool do_hcalin_cluster = false;
0050   bool do_hcalin_eval = false;
0051
0052   bool do_magnet = false;
0053
0054   bool do_hcalout = false;
0055   bool do_hcalout_cell = false;
0056   bool do_hcalout_twr = false;
0057   bool do_hcalout_cluster = false;
0058   bool do_hcalout_eval = false;
0059
0060   bool do_global = false;
0061   bool do_global_fastsim = false;
0062
0063   bool do_jet_reco = false;
0064   bool do_jet_eval = false;
0065
0066   bool do_fwd_jet_reco = false;
0067   bool do_fwd_jet_eval = false;
0068
0069   // fsPHENIX geometry
0070
0071   bool do_FGEM = true;
0072     bool do_FGEM_track =true;
0073     bool do_FGEM_eval = true;
0074
0075   bool do_FEMC = false;
0076   bool do_FEMC_cell = false;
0077   bool do_FEMC_twr = false;
0078   bool do_FEMC_cluster = false;
0079
0080   bool do_FHCAL = false;
0081   bool do_FHCAL_cell = false;
0082   bool do_FHCAL_twr = false;
0083   bool do_FHCAL_cluster = false;
0084
0085   bool do_dst_compress = false;
0086
0087   //Option to convert DST to human command readable TTree for quick poke around the outputs
0088   bool do_DSTReader = true;
0089   //---------------
0090   // Load libraries
0091   //---------------
0092
0093   gSystem->Load("libfun4all.so");
0094   gSystem->Load("libg4detectors.so");
0095   gSystem->Load("libphhepmc.so");
0096   gSystem->Load("libg4testbench.so");
0097   gSystem->Load("libg4hough.so");
0098   gSystem->Load("libg4calo.so");
0099   gSystem->Load("libg4eval.so");
0100
0101   // establish the geometry and reconstruction setup
0102   gROOT->LoadMacro("G4Setup_fsPHENIX.C");
0103   G4Init(do_svtx,do_preshower,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,do_FGEM,do_FEMC,do_FHCAL);
0104
0105   int absorberactive = 0; // set to 1 to make all absorbers active volumes
0106   //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0107   const string magfield = "fieldmap.root"; // fsPHENIX field map by Cesar Luiz da Silva <slash@rcf.rhic.bnl.gov>, sPHENIX + piston
0108   const float magfield_rescale = 1.0; // already adjusted to 1.4T central field
0109
0110   //---------------
0111   // Fun4All server
0112   //---------------
0113
0114   Fun4AllServer *se = Fun4AllServer::instance();
0115 //  se->Verbosity(0); // uncomment for batch production running with minimal output messages
0116   se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running
0117   // just if we set some flags somewhere in this macro
0118   recoConsts *rc = recoConsts::instance();
0119   // By default every random number generator uses
0120   // PHRandomSeed() which reads /dev/urandom to get its seed
0121   // if the RANDOMSEED flag is set its value is taken as seed
0122   // You can either set this to a random value using PHRandomSeed()
0123   // which will make all seeds identical (not sure what the point of
0124   // this would be:
0125   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0126   // or set it to a fixed value so you can debug your code
0127   // rc->set_IntFlag("RANDOMSEED", 12345);
0128
0129   //-----------------
0130   // Event generation
0131   //-----------------
0132
0133   if (readhits)
0134     {
0135       // Get the hits from a file
0136       // The input manager is declared later
0137     }
0138   else if (readhepmc)
0139     {
0140       // this module is needed to read the HepMC records into our G4 sims
0141       // but only if you read HepMC input files
0142       HepMCNodeReader *hr = new HepMCNodeReader();
0143       se->registerSubsystem(hr);
0144     }
0145   else if (runpythia8)
0146     {
0147       gSystem->Load("libPHPythia8.so");
0148
0149       PHPythia8* pythia8 = new PHPythia8();
0150       // see coresoftware/generators/PHPythia8 for example config
0151       pythia8->set_config_file("phpythia8.cfg");
0152       se->registerSubsystem(pythia8);
0153
0154       HepMCNodeReader *hr = new HepMCNodeReader();
0155       se->registerSubsystem(hr);
0156     }
0157   else if (runpythia6)
0158     {
0159       gSystem->Load("libPHPythia6.so");
0160
0161       PHPythia6 *pythia6 = new PHPythia6();
0162       pythia6->set_config_file("phpythia6.cfg");
0163       se->registerSubsystem(pythia6);
0164
0165       HepMCNodeReader *hr = new HepMCNodeReader();
0166       se->registerSubsystem(hr);
0167     }
0168   else
0169     {
0170       // toss low multiplicity dummy events
0171       PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0172       //gen->add_particles("e-",5); // mu+,e+,proton,pi+,Upsilon
0173       //gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
0174       gen->add_particles("mu-",1); // mu-,e-,anti_proton,pi-
0175       //gen->add_particles("pi+",1); // mu-,e-,anti_proton,pi-
0176       if (readhepmc) {
0177     gen->set_reuse_existing_vertex(true);
0178     gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
0179       } else {
0180     gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0181                            PHG4SimpleEventGenerator::Uniform,
0182                            PHG4SimpleEventGenerator::Uniform);
0183     gen->set_vertex_distribution_mean(0.0,0.0,0.0);
0184     gen->set_vertex_distribution_width(0.0,0.0,0.0);
0185       }
0186       gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0187       gen->set_vertex_size_parameters(0.0,0.0);
0188       //gen->set_eta_range(1.4, 3.0);
0189       gen->set_eta_range(3.0, 3.0); //fsPHENIX FWD
0190       gen->set_phi_range(-1.0*TMath::Pi(), 1.0*TMath::Pi());
0191       //gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
0192       gen->set_p_range(0, 30.0);
0193       gen->Embed(1);
0194       gen->Verbosity(0);
0195       se->registerSubsystem(gen);
0196     }
0197
0198   if (!readhits)
0199     {
0200       //---------------------
0201       // Detector description
0202       //---------------------
0203
0204       G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
0205           do_svtx, do_preshower, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe,
0206           do_FGEM, do_FEMC, do_FHCAL,
0207           magfield_rescale);
0208
0209     }
0210
0211   //---------
0212   // BBC Reco
0213   //---------
0214
0215   if (do_bbc)
0216     {
0217       gROOT->LoadMacro("G4_Bbc.C");
0218       BbcInit();
0219       Bbc_Reco();
0220     }
0221
0222   //------------------
0223   // Detector Division
0224   //------------------
0225
0226   if (do_svtx_cell) Svtx_Cells();
0227
0228   if (do_cemc_cell) CEMC_Cells();
0229
0230   if (do_hcalin_cell) HCALInner_Cells();
0231
0232   if (do_hcalout_cell) HCALOuter_Cells();
0233
0234   if (do_FEMC_cell) FEMC_Cells();
0235   if (do_FHCAL_cell) FHCAL_Cells();
0236
0237   //-----------------------------
0238   // CEMC towering and clustering
0239   //-----------------------------
0240
0241   if (do_cemc_twr) CEMC_Towers();
0242   if (do_cemc_cluster) CEMC_Clusters();
0243
0244   //-----------------------------
0245   // HCAL towering and clustering
0246   //-----------------------------
0247
0248   if (do_hcalin_twr) HCALInner_Towers();
0249   if (do_hcalin_cluster) HCALInner_Clusters();
0250
0251   if (do_hcalout_twr) HCALOuter_Towers();
0252   if (do_hcalout_cluster) HCALOuter_Clusters();
0253
0254   if (do_FEMC_twr) FEMC_Towers();
0255   if (do_FEMC_cluster) FEMC_Clusters();
0256
0257   if (do_FHCAL_twr) FHCAL_Towers();
0258   if (do_FHCAL_cluster) FHCAL_Clusters();
0259
0260   if (do_dst_compress) ShowerCompress();
0261
0262   //--------------
0263   // SVTX tracking
0264   //--------------
0265
0266   if (do_svtx_track) Svtx_Reco();
0267
0268   //--------------
0269   // FGEM tracking
0270   //--------------
0271
0272     if(do_FGEM_track) FGEM_FastSim_Reco();
0273
0274   //-----------------
0275   // Global Vertexing
0276   //-----------------
0277
0278   if (do_global)
0279     {
0280       gROOT->LoadMacro("G4_Global.C");
0281       Global_Reco();
0282     }
0283
0284   else if (do_global_fastsim)
0285     {
0286       gROOT->LoadMacro("G4_Global.C");
0287       Global_FastSim();
0288     }
0289
0290   //---------
0291   // Jet reco
0292   //---------
0293
0294   if (do_jet_reco)
0295     {
0296       gROOT->LoadMacro("G4_Jets.C");
0297       Jet_Reco();
0298     }
0299
0300   if (do_fwd_jet_reco)
0301     {
0302       gROOT->LoadMacro("G4_FwdJets.C");
0303       Jet_FwdReco();
0304     }
0305   //----------------------
0306   // Simulation evaluation
0307   //----------------------
0308
0309   if (do_svtx_eval) Svtx_Eval("g4svtx_eval.root");
0310
0311   if (do_cemc_eval) CEMC_Eval("g4cemc_eval.root");
0312
0313   if (do_hcalin_eval) HCALInner_Eval("g4hcalin_eval.root");
0314
0315   if (do_hcalout_eval) HCALOuter_Eval("g4hcalout_eval.root");
0316
0317   if (do_jet_eval) Jet_Eval("g4jet_eval.root");
0318
0319   if (do_fwd_jet_eval) Jet_FwdEval("g4fwdjet_eval.root");
0320
0321   if (do_FGEM_eval) FGEM_FastSim_Eval("g4fwdtrack_fastsim_eval.root");
0322
0323   //--------------
0324   // IO management
0325   //--------------
0326
0327   if (readhits)
0328     {
0329       // Hits file
0330       Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
0331       hitsin->fileopen(inputFile);
0332       se->registerInputManager(hitsin);
0333     }
0334   if (readhepmc)
0335     {
0336       Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0337       se->registerInputManager( in );
0338       se->fileopen( in->Name().c_str(), inputFile );
0339     }
0340   else
0341     {
0342       // for single particle generators we just need something which drives
0343       // the event loop, the Dummy Input Mgr does just that
0344       Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0345       se->registerInputManager( in );
0346     }
0347
0348   if (do_DSTReader)
0349     {
0350       //Convert DST to human command readable TTree for quick poke around the outputs
0351       gROOT->LoadMacro("G4_DSTReader_fsPHENIX.C");
0352
0353       G4DSTreader_fsPHENIX( outputFile, //
0354           /*int*/ absorberactive ,
0355           /*bool*/ do_svtx ,
0356           /*bool*/ do_preshower ,
0357           /*bool*/ do_cemc ,
0358           /*bool*/ do_hcalin ,
0359           /*bool*/ do_magnet ,
0360           /*bool*/ do_hcalout ,
0361           /*bool*/ do_cemc_twr ,
0362           /*bool*/ do_hcalin_twr ,
0363           /*bool*/ do_magnet  ,
0364                 /*bool*/ do_hcalout_twr,
0365                 /*bool*/ do_FGEM,
0366                 /*bool*/ do_FHCAL,
0367                 /*bool*/ do_FHCAL_twr,
0368                 /*bool*/ do_FEMC,
0369                 /*bool*/ do_FEMC_twr
0370           );
0371     }
0372
0373   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0374   if (do_dst_compress) DstCompress(out);
0375   se->registerOutputManager(out);
0376
0377     if (nEvents == 0 && !readhits && !readhepmc)
0378     {
0379         cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0380         cout << "it will run forever, so I just return without running anything" << endl;
0381         return;
0382     }
0383
0384     if (nEvents < 0)
0385     {
0386         PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0387         g4->ApplyCommand("/control/execute vis.mac");
0388         //g4->StartGui();
0389         se->run(1);
0390
0391         se->End();
0392         std::cout << "All done" << std::endl;
0393
0394
0395         std::cout << "==== Useful display commands ==" << std::endl;
0396         cout << "draw axis: " << endl;
0397         cout << " G4Cmd(\"/vis/scene/add/axes 0 0 0 50 cm\")" << endl;
0398         cout << "zoom" << endl;
0399         cout << " G4Cmd(\"/vis/viewer/zoom 1\")" << endl;
0400         cout << "viewpoint:" << endl;
0401         cout << " G4Cmd(\"/vis/viewer/set/viewpointThetaPhi 0 0\")" << endl;
0402         cout << "panTo:" << endl;
0403         cout << " G4Cmd(\"/vis/viewer/panTo 0 0 cm\")" << endl;
0404         cout << "print to eps:" << endl;
0405         cout << " G4Cmd(\"/vis/ogl/printEPS\")" << endl;
0406         cout << "set background color:" << endl;
0407         cout << " G4Cmd(\"/vis/viewer/set/background white\")" << endl;
0408         std::cout << "===============================" << std::endl;
0409     }
0410     else
0411     {
0412
0413         se->run(nEvents);
0414
0415         se->End();
0416
0417         //cin.get();
0418
0419         std::cout << "All done" << std::endl;
0420         delete se;
0421         gSystem->Exit(0);
0422     }
0423
0424 }
0425
0426
0427     void
0428 G4Cmd(const char * cmd)
0429 {
0430     Fun4AllServer *se = Fun4AllServer::instance();
0431     PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0432     g4->ApplyCommand(cmd);
0433 }