macros/leptoquarks_fun4all/Fun4All_G4_EICDetector_LQ.C

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