macros/leptoquarks_fun4all/Fun4All_G4_EICDetector_LQ_reference.C

0001 int Fun4All_G4_EICDetector_LQ_reference(
0002                   string n="1093",
0003                   string ebeam="20",
0004                   string pbeam="250",
0005                   //string inputFile,
0006                   string inputFile="/direct/phenix+u/spjeffas/LQGENEP/TestOut.1093event.root",
0007                   string output="",
0008                   const char * outputFile = "G4EICDetector.root"
0009                )
0010 {
0011   // Set the number of TPC layer
0012   const int n_TPC_layers = 40;  // use 60 for backward compatibility only
0013
0014   //Get parameter variables from parameter file
0015
0016   int nEvents;
0017   stringstream geek(n);
0018   geek>>nEvents;
0019
0020
0021
0022   string directory = "/direct/phenix+u/spjeffas/leptoquark/output/"+output+"/";
0023
0024
0025   //===============
0026   // Input options
0027   //===============
0028   // Either:
0029   // read previously generated g4-hits files, in this case it opens a DST and skips
0030   // the simulations step completely. The G4Setup macro is only loaded to get information
0031   // about the number of layers used for the cell reco code
0032   //
0033   // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
0034   // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0035   const bool readhits = false;
0036   // Or:
0037   // read files in HepMC format (typically output from event generators like hijing or pythia)
0038   const bool readhepmc = false; // read HepMC files
0039   // Or:
0040   // read files in EICTree format generated by eicsmear package
0041   const bool readeictree = true;
0042   // Or:
0043   // Use Pythia 8
0044   const bool runpythia8 = false;
0045   // Or:
0046   // Use Pythia 6
0047   const bool runpythia6 = false;
0048   // Or:
0049   // Use HEPGen
0050   const bool runhepgen = false;
0051   // Or:
0052   // Use Sartre
0053   const bool runsartre = false;
0054
0055
0056
0057   // Besides the above flags. One can further choose to further put in following particles in Geant4 simulation
0058   // Use multi-particle generator (PHG4SimpleEventGenerator), see the code block below to choose particle species and kinematics
0059   const bool particles = false && !readhits;
0060   // or gun/ very simple single particle gun generator
0061   const bool usegun = false && !readhits;
0062   // Throw single Upsilons, may be embedded in Hijing by setting readhepmc flag also  (note, careful to set Z vertex equal to Hijing events)
0063   const bool upsilons = false && !readhits;
0064
0065   //======================
0066   // What to run
0067   //======================
0068
0069   // sPHENIX barrel
0070   bool do_bbc = true;
0071
0072   bool do_pipe = true;
0073
0074   bool do_svtx = true;
0075   bool do_svtx_cell = do_svtx && true;
0076   bool do_svtx_track = do_svtx_cell && true;
0077   bool do_svtx_eval = do_svtx_track && true;
0078
0079   bool do_pstof = false;
0080
0081   bool do_cemc = true;
0082
0083   bool do_cemc_cell = true;
0084   bool do_cemc_twr = true;
0085   bool do_cemc_cluster = true;
0086   bool do_cemc_eval = true;
0087
0088   bool do_hcalin = true;
0089   bool do_hcalin_cell = true;
0090   bool do_hcalin_twr = true;
0091   bool do_hcalin_cluster = true;
0092   bool do_hcalin_eval = true;
0093
0094   bool do_cemc_cell = do_cemc && true;
0095   bool do_cemc_twr = do_cemc_cell && true;
0096   bool do_cemc_cluster = do_cemc_twr && true;
0097   bool do_cemc_eval = do_cemc_cluster && true;
0098
0099   bool do_hcalin = true;
0100   bool do_hcalin_cell = do_hcalin && true;
0101   bool do_hcalin_twr = do_hcalin_cell && true;
0102   bool do_hcalin_cluster = do_hcalin_twr && true;
0103   bool do_hcalin_eval = do_hcalin_cluster && true;
0104
0105
0106   bool do_magnet = true;
0107
0108   bool do_hcalout = true;
0109
0110   bool do_hcalout_cell = true;
0111   bool do_hcalout_twr = true;
0112   bool do_hcalout_cluster = true;
0113   bool do_hcalout_eval = true;
0114
0115   bool do_global = true;
0116   bool do_global_fastsim = false;
0117
0118   bool do_jet_reco = true;
0119   bool do_jet_eval = true;
0120
0121   bool do_fwd_jet_reco = true;
0122   bool do_fwd_jet_eval = false;
0123
0124   bool do_hcalout_cell = do_hcalout && true;
0125   bool do_hcalout_twr = do_hcalout_cell && true;
0126   bool do_hcalout_cluster = do_hcalout_twr && true;
0127   bool do_hcalout_eval = do_hcalout_cluster && true;
0128
0129
0130   // EICDetector geometry - barrel
0131   bool do_DIRC = true;
0132
0133   // EICDetector geometry - 'hadron' direction
0134   bool do_FGEM = true;
0135   bool do_FGEM_track = do_FGEM &&  false;
0136
0137   bool do_RICH = true;
0138   bool do_Aerogel = true;
0139
0140   bool do_FEMC = true;
0141   bool do_FEMC_cell = do_FEMC && true;
0142   bool do_FEMC_twr = do_FEMC_cell && true;
0143   bool do_FEMC_cluster = do_FEMC_twr && true;
0144   bool do_FEMC_eval = do_FEMC_cluster && true;
0145
0146   bool do_FHCAL = true;
0147   bool do_FHCAL_cell = do_FHCAL && true;
0148   bool do_FHCAL_twr = do_FHCAL_cell && true;
0149   bool do_FHCAL_cluster = do_FHCAL_twr && true;
0150   bool do_FHCAL_eval = do_FHCAL_cluster && true;
0151
0152   // EICDetector geometry - 'electron' direction
0153   bool do_EGEM = true;
0154   bool do_EGEM_track = do_EGEM &&  false;
0155
0156   bool do_EEMC = true;
0157   bool do_EEMC_cell = do_EEMC && true;
0158   bool do_EEMC_twr = do_EEMC_cell && true;
0159   bool do_EEMC_cluster = do_EEMC_twr && true;
0160   bool do_EEMC_eval = do_EEMC_cluster && true;
0161
0162   //do leptoquark analysis modules
0163   bool do_lepto_analysis = true;
0164
0165   // Other options
0166   bool do_global = true;
0167   bool do_global_fastsim = false;
0168
0169   bool do_calotrigger = false && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
0170
0171   bool do_jet_reco = true;
0172   bool do_jet_eval = do_jet_reco && true;
0173
0174   bool do_fwd_jet_reco = true;
0175   bool do_fwd_jet_eval = do_fwd_jet_reco && true;
0176
0177   // HI Jet Reco for jet simulations in Au+Au (default is false for
0178   // single particle / p+p simulations, or for Au+Au simulations which
0179   // don't care about jets)
0180   bool do_HIjetreco = false && do_jet_reco && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
0181
0182   // Compress DST files
0183   bool do_dst_compress = false;
0184
0185   //Option to convert DST to human command readable TTree for quick poke around the outputs
0186   bool do_DSTReader = false;
0187
0188   //---------------
0189   // Load libraries
0190   //---------------
0191
0192   gSystem->Load("libfun4all.so");
0193   gSystem->Load("libg4detectors.so");
0194   gSystem->Load("libphhepmc.so");
0195   gSystem->Load("libg4testbench.so");
0196   gSystem->Load("libg4hough.so");
0197   gSystem->Load("libg4calo.so");
0198   gSystem->Load("libg4eval.so");
0199   gSystem->Load("libeicana.so");
0200
0201   // establish the geometry and reconstruction setup
0202   gROOT->LoadMacro("G4Setup_EICDetector.C");
0203   G4Init(do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,do_FGEM,do_EGEM,do_FEMC,do_FHCAL,do_EEMC,do_DIRC,do_RICH,do_Aerogel,n_TPC_layers);
0204
0205   int absorberactive = 0; // set to 1 to make all absorbers active volumes
0206   //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0207   const string magfield = "/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0208   const float magfield_rescale = 1.4/1.5; // scale the map to a 1.4 T field
0209
0210   //---------------
0211   // Fun4All server
0212   //---------------
0213
0214   Fun4AllServer *se = Fun4AllServer::instance();
0215   se->Verbosity(0); // uncomment for batch production running with minimal output messages
0216   // se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running
0217
0218   // just if we set some flags somewhere in this macro
0219   recoConsts *rc = recoConsts::instance();
0220   // By default every random number generator uses
0221   // PHRandomSeed() which reads /dev/urandom to get its seed
0222   // if the RANDOMSEED flag is set its value is taken as seed
0223   // You can either set this to a random value using PHRandomSeed()
0224   // which will make all seeds identical (not sure what the point of
0225   // this would be:
0226   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0227   // or set it to a fixed value so you can debug your code
0228   // rc->set_IntFlag("RANDOMSEED", 12345);
0229
0230   //-----------------
0231   // Event generation
0232   //-----------------
0233
0234   if (readhits)
0235     {
0236       // Get the hits from a file
0237       // The input manager is declared later
0238     }
0239   else if (readhepmc)
0240     {
0241       // this module is needed to read the HepMC records into our G4 sims
0242       // but only if you read HepMC input files
0243       HepMCNodeReader *hr = new HepMCNodeReader();
0244       se->registerSubsystem(hr);
0245     }
0246   else if (readeictree)
0247     {
0248       // this module is needed to read the EICTree style records into our G4 sims
0249       ReadEICFiles *eicr = new ReadEICFiles();
0250       eicr->OpenInputFile(inputFile);
0251
0252       se->registerSubsystem(eicr);
0253     }
0254   else if (runpythia8)
0255     {
0256       gSystem->Load("libPHPythia8.so");
0257
0258       PHPythia8* pythia8 = new PHPythia8();
0259       // see coresoftware/generators/PHPythia8 for example config
0260       pythia8->set_config_file("/direct/phenix+u/spjeffas/coresoftware/generators/PHPythia8/phpythia8.cfg");
0261       se->registerSubsystem(pythia8);
0262
0263       HepMCNodeReader *hr = new HepMCNodeReader();
0264       se->registerSubsystem(hr);
0265     }
0266   else if (runpythia6)
0267     {
0268       gSystem->Load("libPHPythia6.so");
0269
0270       PHPythia6 *pythia6 = new PHPythia6();
0271       // see coresoftware/generators/PHPythia6 for example config
0272       pythia6->set_config_file("/direct/phenix+u/spjeffas/coresoftware/generators/PHPythia6/phpythia6_ep.cfg");
0273       se->registerSubsystem(pythia6);
0274
0275       HepMCNodeReader *hr = new HepMCNodeReader();
0276       se->registerSubsystem(hr);
0277     }
0278   else if (runhepgen)
0279     {
0280       gSystem->Load("libsHEPGen.so");
0281
0282       sHEPGen *hepgen = new sHEPGen();
0283       // see HEPGen source directory/share/vggdata for required .dat files
0284       // see HEPGen source directory/share/datacards for required datacard files
0285       hepgen->set_datacard_file("hepgen_dvcs.data");
0286       hepgen->set_momentum_electron(-20);
0287       hepgen->set_momentum_hadron(250);
0288       se->registerSubsystem(hepgen);
0289
0290       HepMCNodeReader *hr = new HepMCNodeReader();
0291       se->registerSubsystem(hr);
0292     }
0293   else if (runsartre)
0294     {
0295       // see coresoftware/generators/PHSartre/README for setup instructions
0296       // before running:
0297       // setenv SARTRE_DIR /opt/sphenix/core/sartre-1.20_root-5.34.36
0298       gSystem->Load("libPHSartre.so");
0299
0300       PHSartre* mysartre = new PHSartre();
0301       // see coresoftware/generators/PHSartre for example config
0302       mysartre->set_config_file("sartre.cfg");
0303
0304       // particle trigger to enhance forward J/Psi -> ee
0305       PHSartreParticleTrigger* pTrig = new PHSartreParticleTrigger("MySartreTrigger");
0306       pTrig->AddParticles(-11);
0307       //pTrig->SetEtaHighLow(4.0,1.4);
0308       pTrig->SetEtaHighLow(1.0,-1.1);  // central arm
0309       pTrig->PrintConfig();
0310       mysartre->register_trigger((PHSartreGenTrigger *)pTrig);
0311       se->registerSubsystem(mysartre);
0312
0313       HepMCNodeReader *hr = new HepMCNodeReader();
0314       se->registerSubsystem(hr);
0315     }
0316
0317   // If "readhepMC" is also set, the particles will be embedded in Hijing events
0318   if(particles)
0319     {
0320       // toss low multiplicity dummy events
0321       PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0322
0323       //gen->add_particles("e-",5); // mu+,e+,proton,pi+,Upsilon
0324       //gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
0325       gen->add_particles("tau-",1); // mu-,e-,anti_proton,pi-
0326       if (readhepmc) {
0327         gen->set_reuse_existing_vertex(true);
0328         gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
0329       } else {
0330         gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0331                                               PHG4SimpleEventGenerator::Uniform,
0332                                               PHG4SimpleEventGenerator::Uniform);
0333         gen->set_vertex_distribution_mean(0.0,0.0,0.0);
0334         gen->set_vertex_distribution_width(0.0,0.0,5.0);
0335       }
0336       gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0337       gen->set_vertex_size_parameters(0.0,0.0);
0338       gen->set_eta_range(0.1, 0.1);
0339       //gen->set_eta_range(3.0, 3.0); //EICDetector FWD
0340       gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
0341       //gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
0342       gen->set_p_range(30.0, 30.0);
0343
0344       //gen->add_particles("pi-",1); // mu+,e+,proton,pi+,Upsilon
0345       //gen->add_particles("pi+",100); // 100 pion option
0346       if (readhepmc)
0347         {
0348           gen->set_reuse_existing_vertex(true);
0349           gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0350         }
0351       else
0352         {
0353           gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0354                                                 PHG4SimpleEventGenerator::Uniform,
0355                                                 PHG4SimpleEventGenerator::Uniform);
0356           gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
0357           gen->set_vertex_distribution_width(0.0, 0.0, 0.0);
0358         }
0359       gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0360       gen->set_vertex_size_parameters(0.0, 0.0);
0361       gen->set_eta_range(-1.0, 1.0);
0362       gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0363       //gen->set_pt_range(0.1, 50.0);
0364       gen->set_pt_range(0.1, 20.0);
0365
0366       gen->Embed(1);
0367       gen->Verbosity(0);
0368
0369       se->registerSubsystem(gen);
0370     }
0371   if (usegun)
0372     {
0373       // PHG4ParticleGun *gun = new PHG4ParticleGun();
0374       // gun->set_name("anti_proton");
0375       // gun->set_name("geantino");
0376       // gun->set_vtx(0, 0, 0);
0377       // gun->set_mom(10, 0, 0.01);
0378       // gun->AddParticle("geantino",1.7776,-0.4335,0.);
0379       // gun->AddParticle("geantino",1.7709,-0.4598,0.);
0380       // gun->AddParticle("geantino",2.5621,0.60964,0.);
0381       // gun->AddParticle("geantino",1.8121,0.253,0.);
0382       // se->registerSubsystem(gun);
0383       PHG4ParticleGenerator *pgen = new PHG4ParticleGenerator();
0384       pgen->set_name("e-");
0385       pgen->set_z_range(0,0);
0386       pgen->set_eta_range(0.01,0.01);
0387       pgen->set_mom_range(10,10);
0388       pgen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0389       se->registerSubsystem(pgen);
0390     }
0391
0392   // If "readhepMC" is also set, the Upsilons will be embedded in Hijing events, if 'particles" is set, the Upsilons will be embedded in whatever particles are thrown
0393   if(upsilons)
0394     {
0395       // run upsilons for momentum, dca performance, alone or embedded in Hijing
0396
0397       PHG4ParticleGeneratorVectorMeson *vgen = new PHG4ParticleGeneratorVectorMeson();
0398       vgen->add_decay_particles("e+","e-",0); // i = decay id
0399       // event vertex
0400       if (readhepmc || particles)
0401         {
0402           vgen->set_reuse_existing_vertex(true);
0403         }
0404       else
0405         {
0406           vgen->set_vtx_zrange(-10.0, +10.0);
0407         }
0408
0409       // Note: this rapidity range completely fills the acceptance of eta = +/- 1 unit
0410       vgen->set_rapidity_range(-1.0, +1.0);
0411       vgen->set_pt_range(0.0, 10.0);
0412
0413       int istate = 1;
0414
0415       if(istate == 1)
0416         {
0417           // Upsilon(1S)
0418           vgen->set_mass(9.46);
0419           vgen->set_width(54.02e-6);
0420         }
0421       else if (istate == 2)
0422         {
0423           // Upsilon(2S)
0424           vgen->set_mass(10.0233);
0425           vgen->set_width(31.98e-6);
0426         }
0427       else
0428         {
0429           // Upsilon(3S)
0430           vgen->set_mass(10.3552);
0431           vgen->set_width(20.32e-6);
0432         }
0433
0434       vgen->Verbosity(0);
0435       vgen->Embed(2);
0436       se->registerSubsystem(vgen);
0437
0438       cout << "Upsilon generator for istate = " << istate << " created and registered "  << endl;
0439     }
0440
0441   if (!readhits)
0442     {
0443       //---------------------
0444       // Detector description
0445       //---------------------
0446
0447       G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
0448               do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,
0449               do_FGEM,do_EGEM,do_FEMC,do_FHCAL,do_EEMC,do_DIRC,do_RICH,do_Aerogel,
0450               magfield_rescale);
0451
0452     }
0453
0454   //---------
0455   // BBC Reco
0456   //---------
0457
0458   if (do_bbc)
0459     {
0460       gROOT->LoadMacro("G4_Bbc.C");
0461       BbcInit();
0462       Bbc_Reco();
0463     }
0464
0465   //------------------
0466   // Detector Division
0467   //------------------
0468
0469   if (do_svtx_cell) Svtx_Cells();
0470
0471   if (do_cemc_cell) CEMC_Cells();
0472
0473   if (do_hcalin_cell) HCALInner_Cells();
0474
0475   if (do_hcalout_cell) HCALOuter_Cells();
0476
0477   if (do_FEMC_cell) FEMC_Cells();
0478
0479   if (do_FHCAL_cell) FHCAL_Cells();
0480
0481   if (do_EEMC_cell) EEMC_Cells();
0482
0483   //-----------------------------
0484   // CEMC towering and clustering
0485   //-----------------------------
0486
0487   if (do_cemc_twr) CEMC_Towers();
0488   if (do_cemc_cluster) CEMC_Clusters();
0489
0490   //-----------------------------
0491   // HCAL towering and clustering
0492   //-----------------------------
0493
0494   if (do_hcalin_twr) HCALInner_Towers();
0495   if (do_hcalin_cluster) HCALInner_Clusters();
0496
0497   if (do_hcalout_twr) HCALOuter_Towers();
0498   if (do_hcalout_cluster) HCALOuter_Clusters();
0499
0500   //-----------------------------
0501   // e, h direction Calorimeter  towering and clustering
0502   //-----------------------------
0503
0504   if (do_FEMC_twr) FEMC_Towers();
0505   if (do_FEMC_cluster) FEMC_Clusters();
0506
0507   if (do_FHCAL_twr) FHCAL_Towers();
0508   if (do_FHCAL_cluster) FHCAL_Clusters();
0509
0510   if (do_EEMC_twr) EEMC_Towers();
0511   if (do_EEMC_cluster) EEMC_Clusters();
0512
0513   if (do_dst_compress) ShowerCompress();
0514
0515   //--------------
0516   // SVTX tracking
0517   //--------------
0518
0519   if (do_svtx_track) Svtx_Reco();
0520
0521   //--------------
0522   // FGEM tracking
0523   //--------------
0524
0525   if(do_FGEM_track) FGEM_FastSim_Reco();
0526
0527   //--------------
0528   // EGEM tracking
0529   //--------------
0530
0531   if(do_EGEM_track) EGEM_FastSim_Reco();
0532
0533   //-----------------
0534   // Global Vertexing
0535   //-----------------
0536
0537   if (do_global)
0538     {
0539       gROOT->LoadMacro("G4_Global.C");
0540       Global_Reco();
0541     }
0542
0543   else if (do_global_fastsim)
0544     {
0545       gROOT->LoadMacro("G4_Global.C");
0546       Global_FastSim();
0547     }
0548
0549   //-----------------
0550   // Calo Trigger Simulation
0551   //-----------------
0552
0553   if (do_calotrigger)
0554     {
0555       gROOT->LoadMacro("G4_CaloTrigger.C");
0556       CaloTrigger_Sim();
0557     }
0558
0559   //---------
0560   // Jet reco
0561   //---------
0562
0563   if (do_jet_reco)
0564     {
0565       gROOT->LoadMacro("G4_Jets.C");
0566       Jet_Reco();
0567     }
0568
0569   if (do_HIjetreco) {
0570     gROOT->LoadMacro("G4_HIJetReco.C");
0571     HIJetReco();
0572   }
0573
0574   if (do_fwd_jet_reco)
0575     {
0576       gROOT->LoadMacro("G4_FwdJets.C");
0577       Jet_FwdReco();
0578     }
0579
0580   //----------------------
0581   // Simulation evaluation
0582   //----------------------
0583
0584
0585   if (do_svtx_eval) Svtx_Eval(directory+"g4svtx_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
0586
0587   if (do_cemc_eval) CEMC_Eval(directory+"g4cemc_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
0588
0589   if (do_hcalin_eval) HCALInner_Eval(directory+"g4hcalin_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
0590
0591   if (do_hcalout_eval) HCALOuter_Eval(directory+"g4hcalout_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
0592
0593   if (do_jet_eval) Jet_Eval(directory+"g4jet_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
0594
0595   if (do_fwd_jet_eval) Jet_FwdEval(directory+"g4fwdjet_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
0596
0597   if(do_lepto_analysis){
0598     gROOT->LoadMacro("G4_Lepto.C");
0599     G4_Lepto(directory+"LeptoAna_p"+pbeam+"_e"+ebeam+"_"+n+"events");
0600   }
0601
0602
0603
0604
0605
0606   //--------------
0607   // IO management
0608   //--------------
0609
0610   if (readhits)
0611     {
0612       // Hits file
0613       Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
0614       hitsin->fileopen(inputFile);
0615       se->registerInputManager(hitsin);
0616     }
0617   if (readhepmc)
0618     {
0619       Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0620       se->registerInputManager( in );
0621       se->fileopen( in->Name().c_str(), inputFile );
0622     }
0623   else
0624     {
0625       // for single particle generators we just need something which drives
0626       // the event loop, the Dummy Input Mgr does just that
0627       Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0628       se->registerInputManager( in );
0629     }
0630
0631   if (do_DSTReader)
0632     {
0633       //Convert DST to human command readable TTree for quick poke around the outputs
0634       gROOT->LoadMacro("G4_DSTReader_EICDetector.C");
0635
0636       G4DSTreader_EICDetector( outputFile, //
0637                                /*int*/ absorberactive ,
0638                                /*bool*/ do_svtx ,
0639                                /*bool*/ do_cemc ,
0640                                /*bool*/ do_hcalin ,
0641                                /*bool*/ do_magnet ,
0642                                /*bool*/ do_hcalout ,
0643                                /*bool*/ do_cemc_twr ,
0644                                /*bool*/ do_hcalin_twr ,
0645                                /*bool*/ do_magnet  ,
0646                                /*bool*/ do_hcalout_twr,
0647                                /*bool*/ do_FGEM,
0648                                /*bool*/ do_EGEM,
0649                                /*bool*/ do_FHCAL,
0650                                /*bool*/ do_FHCAL_twr,
0651                                /*bool*/ do_FEMC,
0652                                /*bool*/ do_FEMC_twr,
0653                                /*bool*/ do_EEMC,
0654                                /*bool*/ do_EEMC_twr
0655                                );
0656     }
0657
0658   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0659   if (do_dst_compress) DstCompress(out);
0660   se->registerOutputManager(out);
0661
0662   //-----------------
0663   // Event processing
0664   //-----------------
0665   if (nEvents < 0)
0666     {
0667       return;
0668     }
0669   // if we run the particle generator and use 0 it'll run forever
0670   if (nEvents == 0 && !readhits && !readhepmc)
0671     {
0672       cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0673       cout << "it will run forever, so I just return without running anything" << endl;
0674       return;
0675     }
0676
0677   se->run(nEvents);
0678
0679   //-----
0680   // Exit
0681   //-----
0682
0683   se->End();
0684   std::cout << "All done" << std::endl;
0685   delete se;
0686   gSystem->Exit(0);
0687 }