AnaTutorial/macro/Fun4All_AnaTutorial_sPHENIX.C

0001 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0002 #define MACRO_FUN4ALLG4SPHENIX_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <DisplayOn.C>
0007 #include <G4Setup_sPHENIX.C>
0008 #include <G4_Mbd.C>
0009 #include <G4_CaloTrigger.C>
0010 #include <G4_Centrality.C>
0011 #include <G4_DSTReader.C>
0012 #include <G4_Global.C>
0013 #include <G4_HIJetReco.C>
0014 #include <G4_Input.C>
0015 #include <G4_Jets.C>
0016 #include <G4_KFParticle.C>
0017 #include <G4_ParticleFlow.C>
0018 #include <G4_Production.C>
0019 #include <G4_TopoClusterReco.C>
0020
0021 #include <Trkr_RecoInit.C>
0022 #include <Trkr_Clustering.C>
0023 #include <Trkr_LaserClustering.C>
0024 #include <Trkr_Reco.C>
0025 #include <Trkr_Eval.C>
0026 #include <Trkr_QA.C>
0027
0028 #include <Trkr_Diagnostics.C>
0029 #include <G4_User.C>
0030 #include <QA.C>
0031
0032 #include <anatutorial/AnaTutorial.h>
0033
0034 #include <ffamodules/FlagHandler.h>
0035 #include <ffamodules/HeadReco.h>
0036 #include <ffamodules/SyncReco.h>
0037 #include <ffamodules/CDBInterface.h>
0038
0039 #include <fun4all/Fun4AllDstOutputManager.h>
0040 #include <fun4all/Fun4AllOutputManager.h>
0041 #include <fun4all/Fun4AllServer.h>
0042
0043 #include <phool/PHRandomSeed.h>
0044 #include <phool/recoConsts.h>
0045
0046 R__LOAD_LIBRARY(libfun4all.so)
0047 R__LOAD_LIBRARY(libffamodules.so)
0048
0049 R__LOAD_LIBRARY(libanatutorial.so)
0050 // For HepMC Hijing
0051 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0052
0053 int Fun4All_AnaTutorial_sPHENIX(
0054     const int nEvents = 1,
0055     const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
0056     const string &outputFile = "G4sPHENIX.root",
0057     const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
0058     const int skip = 0,
0059     const string &outdir = ".")
0060 {
0061   Fun4AllServer *se = Fun4AllServer::instance();
0062   se->Verbosity(0);
0063
0064   //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0065   PHRandomSeed::Verbosity(1);
0066
0067   // just if we set some flags somewhere in this macro
0068   recoConsts *rc = recoConsts::instance();
0069   // By default every random number generator uses
0070   // PHRandomSeed() which reads /dev/urandom to get its seed
0071   // if the RANDOMSEED flag is set its value is taken as seed
0072   // You can either set this to a random value using PHRandomSeed()
0073   // which will make all seeds identical (not sure what the point of
0074   // this would be:
0075   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0076   // or set it to a fixed value so you can debug your code
0077   //  rc->set_IntFlag("RANDOMSEED", 12345);
0078
0079
0080   //===============
0081   // Input options
0082   //===============
0083   // verbosity setting (applies to all input managers)
0084   Input::VERBOSITY = 0;
0085   // First enable the input generators
0086   // Either:
0087   // read previously generated g4-hits files, in this case it opens a DST and skips
0088   // the simulations step completely. The G4Setup macro is only loaded to get information
0089   // about the number of layers used for the cell reco code
0090   //  Input::READHITS = true;
0091   INPUTREADHITS::filename[0] = inputFile;
0092   // if you use a filelist
0093   // INPUTREADHITS::listfile[0] = inputFile;
0094   // Or:
0095   // Use particle generator
0096   // And
0097   // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
0098   // In case embedding into a production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
0099   // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0100   //  Input::EMBED = true;
0101   INPUTEMBED::filename[0] = embed_input_file;
0102   // if you use a filelist
0103   //INPUTEMBED::listfile[0] = embed_input_file;
0104
0105   Input::SIMPLE = true;
0106   // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
0107   // Input::SIMPLE_VERBOSITY = 1;
0108
0109   // Enable this is emulating the nominal pp/pA/AA collision vertex distribution
0110   // Input::BEAM_CONFIGURATION = Input::AA_COLLISION; // Input::AA_COLLISION (default), Input::pA_COLLISION, Input::pp_COLLISION
0111
0112   //  Input::PYTHIA6 = true;
0113
0114   // Input::PYTHIA8 = true;
0115
0116   //  Input::GUN = true;
0117   //  Input::GUN_NUMBER = 3; // if you need 3 of them
0118   // Input::GUN_VERBOSITY = 1;
0119
0120   // Input::COSMIC = true;
0121
0122   //D0 generator
0123   //Input::DZERO = false;
0124   //Input::DZERO_VERBOSITY = 0;
0125   //Lambda_c generator //Not ready yet
0126   //Input::LAMBDAC = false;
0127   //Input::LAMBDAC_VERBOSITY = 0;
0128   // Upsilon generator
0129   //Input::UPSILON = true;
0130   //Input::UPSILON_NUMBER = 3; // if you need 3 of them
0131   //Input::UPSILON_VERBOSITY = 0;
0132
0133   //  Input::HEPMC = true;
0134   INPUTHEPMC::filename = inputFile;
0135   //-----------------
0136   // Hijing options (symmetrize hijing, add flow, add fermi motion)
0137   //-----------------
0138   //  INPUTHEPMC::HIJINGFLIP = true;
0139   //  INPUTHEPMC::FLOW = true;
0140   //  INPUTHEPMC::FLOW_VERBOSITY = 3;
0141   //  INPUTHEPMC::FERMIMOTION = true;
0142
0143
0144   // Event pile up simulation with collision rate in Hz MB collisions.
0145   //Input::PILEUPRATE = 50e3; // 50 kHz for AuAu
0146   //Input::PILEUPRATE = 3e6; // 3MHz for pp
0147
0148   // Enable this is emulating the nominal pp/pA/AA collision vertex distribution
0149   // for HepMC records (hijing, pythia8)
0150   //  Input::BEAM_CONFIGURATION = Input::AA_COLLISION; // for 2023 sims we want the AA geometry for no pileup sims
0151   //  Input::BEAM_CONFIGURATION = Input::pp_COLLISION; // for 2024 sims we want the pp geometry for no pileup sims
0152   //  Input::BEAM_CONFIGURATION = Input::pA_COLLISION; // for pAu sims we want the pA geometry for no pileup sims
0153
0154   //-----------------
0155   // Initialize the selected Input/Event generation
0156   //-----------------
0157   // This creates the input generator(s)
0158   InputInit();
0159
0160   //--------------
0161   // Set generator specific options
0162   //--------------
0163   // can only be set after InputInit() is called
0164
0165   // Simple Input generator:
0166   // if you run more than one of these Input::SIMPLE_NUMBER > 1
0167   // add the settings for other with [1], next with [2]...
0168   if (Input::SIMPLE)
0169   {
0170     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 5);
0171     if (Input::HEPMC || Input::EMBED)
0172     {
0173       INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
0174       INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0175     }
0176     else
0177     {
0178       INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
0179                                                                                 PHG4SimpleEventGenerator::Gaus,
0180                                                                                 PHG4SimpleEventGenerator::Gaus);
0181       INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
0182       INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0.01, 0.01, 5.);
0183     }
0184     INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1, 1);
0185     INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
0186     INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.1, 20.);
0187   }
0188   // Upsilons
0189   // if you run more than one of these Input::UPSILON_NUMBER > 1
0190   // add the settings for other with [1], next with [2]...
0191   if (Input::UPSILON)
0192   {
0193     INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("e", 0);
0194     INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
0195     INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
0196     // Y species - select only one, last one wins
0197     INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
0198     if (Input::HEPMC || Input::EMBED)
0199     {
0200       INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
0201       INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0202     }
0203   }
0204   // particle gun
0205   // if you run more than one of these Input::GUN_NUMBER > 1
0206   // add the settings for other with [1], next with [2]...
0207   if (Input::GUN)
0208   {
0209     INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
0210     INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
0211   }
0212
0213   // pythia6
0214   if (Input::PYTHIA6)
0215   {
0216     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0217     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia6);
0218   }
0219   // pythia8
0220   if (Input::PYTHIA8)
0221   {
0222     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0223     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
0224   }
0225
0226   //--------------
0227   // Set Input Manager specific options
0228   //--------------
0229   // can only be set after InputInit() is called
0230
0231   if (Input::HEPMC)
0232   {
0233     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0234     Input::ApplysPHENIXBeamParameter(INPUTMANAGER::HepMCInputManager);
0235
0236     // optional overriding beam parameters
0237     //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 8, 0);  //optional collision smear in space, time
0238     //    INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
0239     // //optional choice of vertex distribution function in space, time
0240     //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
0241     //! embedding ID for the event
0242     //! positive ID is the embedded event of interest, e.g. jetty event from pythia
0243     //! negative IDs are backgrounds, .e.g out of time pile up collisions
0244     //! Usually, ID = 0 means the primary Au+Au collision background
0245     //INPUTMANAGER::HepMCInputManager->set_embedding_id(Input::EmbedID);
0246     if (Input::PILEUPRATE > 0)
0247     {
0248       // Copy vertex settings from foreground hepmc input
0249       INPUTMANAGER::HepMCPileupInputManager->CopyHelperSettings(INPUTMANAGER::HepMCInputManager);
0250       // and then modify the ones you want to be different
0251       // INPUTMANAGER::HepMCPileupInputManager->set_vertex_distribution_width(100e-4,100e-4,8,0);
0252     }
0253   }
0254   if (Input::PILEUPRATE > 0)
0255   {
0256     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0257     Input::ApplysPHENIXBeamParameter(INPUTMANAGER::HepMCPileupInputManager);
0258   }
0259   // register all input generators with Fun4All
0260   InputRegister();
0261
0262   if (! Input::READHITS)
0263   {
0264     rc->set_IntFlag("RUNNUMBER",1);
0265
0266     SyncReco *sync = new SyncReco();
0267     se->registerSubsystem(sync);
0268
0269     HeadReco *head = new HeadReco();
0270     se->registerSubsystem(head);
0271   }
0272 // Flag Handler is always needed to read flags from input (if used)
0273 // and update our rc flags with them. At the end it saves all flags
0274 // again on the DST in the Flags node under the RUN node
0275   FlagHandler *flag = new FlagHandler();
0276   se->registerSubsystem(flag);
0277
0278   // set up production relatedstuff
0279   //   Enable::PRODUCTION = true;
0280
0281   //======================
0282   // Write the DST
0283   //======================
0284
0285   //Enable::DSTOUT = true;
0286   Enable::DSTOUT_COMPRESS = false;
0287   DstOut::OutputDir = outdir;
0288   DstOut::OutputFile = outputFile;
0289
0290   //Option to convert DST to human command readable TTree for quick poke around the outputs
0291   //  Enable::DSTREADER = true;
0292
0293   // turn the display on (default off)
0294    //Enable::DISPLAY = true;
0295
0296   //======================
0297   // What to run
0298   //======================
0299
0300   // QA, main switch
0301   Enable::QA = true;
0302
0303   // Global options (enabled for all enables subsystems - if implemented)
0304   //  Enable::ABSORBER = true;
0305   //  Enable::OVERLAPCHECK = true;
0306   //  Enable::VERBOSITY = 1;
0307
0308   // Enable::MBD = true;
0309   // Enable::MBD_SUPPORT = true; // save hist in MBD/BBC support structure
0310   // Enable::MBDRECO = Enable::MBD && true;
0311   Enable::MBDFAKE = true;  // Smeared vtx and t0, use if you don't want real MBD/BBC in simulation
0312
0313   Enable::PIPE = true;
0314   Enable::PIPE_ABSORBER = true;
0315
0316   // central tracking
0317   Enable::MVTX = true;
0318   Enable::MVTX_CELL = Enable::MVTX && true;
0319   Enable::MVTX_CLUSTER = Enable::MVTX_CELL && true;
0320   Enable::MVTX_QA = Enable::MVTX_CLUSTER && Enable::QA && true;
0321
0322   Enable::INTT = true;
0323 //  Enable::INTT_ABSORBER = true; // enables layerwise support structure readout
0324 //  Enable::INTT_SUPPORT = true; // enable global support structure readout
0325   Enable::INTT_CELL = Enable::INTT && true;
0326   Enable::INTT_CLUSTER = Enable::INTT_CELL && true;
0327   Enable::INTT_QA = Enable::INTT_CLUSTER && Enable::QA && true;
0328
0329   Enable::TPC = true;
0330   Enable::TPC_ABSORBER = true;
0331   Enable::TPC_CELL = Enable::TPC && true;
0332   Enable::TPC_CLUSTER = Enable::TPC_CELL && true;
0333   Enable::TPC_QA = Enable::TPC_CLUSTER && Enable::QA && true;
0334
0335   Enable::MICROMEGAS = true;
0336   Enable::MICROMEGAS_CELL = Enable::MICROMEGAS && true;
0337   Enable::MICROMEGAS_CLUSTER = Enable::MICROMEGAS_CELL && true;
0338   Enable::MICROMEGAS_QA = Enable::MICROMEGAS_CLUSTER && Enable::QA && true;
0339
0340   Enable::TRACKING_TRACK = (Enable::MICROMEGAS_CLUSTER && Enable::TPC_CLUSTER && Enable::INTT_CLUSTER && Enable::MVTX_CLUSTER) && true;
0341   Enable::GLOBAL_RECO = (Enable::MBDFAKE || Enable::MBDRECO || Enable::TRACKING_TRACK) && true;
0342   Enable::TRACKING_EVAL = Enable::TRACKING_TRACK && Enable::GLOBAL_RECO && true;
0343   Enable::TRACKING_QA = Enable::TRACKING_TRACK && Enable::QA && true;
0344
0345   // only do track matching if TRACKINGTRACK is also used
0346   Enable::TRACK_MATCHING = Enable::TRACKING_TRACK && false;
0347   Enable::TRACK_MATCHING_TREE = Enable::TRACK_MATCHING && false;
0348   Enable::TRACK_MATCHING_TREE_CLUSTERS = Enable::TRACK_MATCHING_TREE && false;
0349
0350   //Additional tracking tools
0351   //Enable::TRACKING_DIAGNOSTICS = Enable::TRACKING_TRACK && true;
0352   //G4TRACKING::filter_conversion_electrons = true;
0353   // G4TRACKING::use_alignment = true;
0354
0355   // enable pp mode and set extended readout time
0356   // TRACKING::pp_mode = true;
0357   // TRACKING::pp_extended_readout_time = 20000;
0358
0359   // set flags to simulate and correct TPC distortions, specify distortion and correction files
0360   //G4TPC::ENABLE_STATIC_DISTORTIONS = true;
0361   //G4TPC::static_distortion_filename = std::string("/sphenix/user/rcorliss/distortion_maps/2023.02/Summary_hist_mdc2_UseFieldMaps_AA_event_0_bX180961051_0.distortion_map.hist.root");
0362   //G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0363   //G4TPC::static_correction_filename = std::string("/sphenix/user/rcorliss/distortion_maps/2023.02/Summary_hist_mdc2_UseFieldMaps_AA_smoothed_average.correction_map.hist.root");
0364   //G4TPC::ENABLE_AVERAGE_CORRECTIONS = false;
0365
0366   //  cemc electronics + thin layer of W-epoxy to get albedo from cemc
0367   //  into the tracking, cannot run together with CEMC
0368   //  Enable::CEMCALBEDO = true;
0369
0370   Enable::CEMC = true;
0371   Enable::CEMC_ABSORBER = true;
0372   Enable::CEMC_CELL = Enable::CEMC && true;
0373   Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
0374   Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
0375   Enable::CEMC_EVAL = Enable::CEMC_G4Hit && Enable::CEMC_CLUSTER && true;
0376   Enable::CEMC_QA = Enable::CEMC_CLUSTER && Enable::QA && true;
0377
0378   Enable::HCALIN = true;
0379   Enable::HCALIN_ABSORBER = true;
0380   Enable::HCALIN_CELL = Enable::HCALIN && true;
0381   Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
0382   Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
0383   Enable::HCALIN_EVAL = Enable::HCALIN_G4Hit && Enable::HCALIN_CLUSTER && true;
0384   Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && true;
0385
0386   Enable::MAGNET = true;
0387   Enable::MAGNET_ABSORBER = true;
0388
0389   Enable::HCALOUT = true;
0390   Enable::HCALOUT_ABSORBER = true;
0391   Enable::HCALOUT_CELL = Enable::HCALOUT && true;
0392   Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
0393   Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
0394   Enable::HCALOUT_EVAL = Enable::HCALOUT_G4Hit && Enable::HCALOUT_CLUSTER && true;
0395   Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && true;
0396
0397   Enable::EPD = true;
0398   Enable::EPD_TILE = Enable::EPD && true;
0399
0400   Enable::BEAMLINE = true;
0401   //  Enable::BEAMLINE_ABSORBER = true;  // makes the beam line magnets sensitive volumes
0402   //  Enable::BEAMLINE_BLACKHOLE = true; // turns the beamline magnets into black holes
0403   Enable::ZDC = true;
0404   //  Enable::ZDC_ABSORBER = true;
0405   //  Enable::ZDC_SUPPORT = true;
0406   Enable::ZDC_TOWER = Enable::ZDC && true;
0407   Enable::ZDC_EVAL = Enable::ZDC_TOWER && true;
0408
0409   //! forward flux return plug door. Out of acceptance and off by default.
0410   //Enable::PLUGDOOR = true;
0411   Enable::PLUGDOOR_ABSORBER = true;
0412
0413  //Enable::GLOBAL_FASTSIM = true;
0414
0415   //Enable::KFPARTICLE = true;
0416   //Enable::KFPARTICLE_VERBOSITY = 1;
0417   //Enable::KFPARTICLE_TRUTH_MATCH = true;
0418   //Enable::KFPARTICLE_SAVE_NTUPLE = true;
0419
0420   Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0421
0422   Enable::JETS = (Enable::GLOBAL_RECO || Enable::GLOBAL_FASTSIM) && true;
0423   Enable::JETS_EVAL = Enable::JETS && true;
0424   Enable::JETS_QA = Enable::JETS && Enable::QA && true;
0425
0426   // HI Jet Reco for p+Au / Au+Au collisions (default is false for
0427   // single particle / p+p-only simulations, or for p+Au / Au+Au
0428   // simulations which don't particularly care about jets)
0429   Enable::HIJETS = Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0430
0431   // 3-D topoCluster reconstruction, potentially in all calorimeter layers
0432   Enable::TOPOCLUSTER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0433   // particle flow jet reconstruction - needs topoClusters!
0434   Enable::PARTICLEFLOW = Enable::TOPOCLUSTER && true;
0435   // centrality reconstruction
0436   Enable::CENTRALITY = true;
0437
0438   // new settings using Enable namespace in GlobalVariables.C
0439   Enable::BLACKHOLE = true;
0440   //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
0441   //Enable::BLACKHOLE_FORWARD_SAVEHITS = false; // disable forward/backward hits
0442   //BlackHoleGeometry::visible = true;
0443
0444   // run user provided code (from local G4_User.C)
0445   //Enable::USER = true;
0446
0447   //===============
0448   // conditions DB flags
0449   //===============
0450   Enable::CDB = true;
0451   // global tag
0452   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0453   // 64 bit timestamp
0454   rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
0455   //---------------
0456   // World Settings
0457   //---------------
0458   //  G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
0459   //  G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
0460
0461   //---------------
0462   // Magnet Settings
0463   //---------------
0464
0465   //  G4MAGNET::magfield =  string(getenv("CALIBRATIONROOT"))+ string("/Field/Map/sphenix3dbigmapxyz.root");  // default map from the calibration database
0466   //  G4MAGNET::magfield = "1.5"; // alternatively to specify a constant magnetic field, give a float number, which will be translated to solenoidal field in T, if string use as fieldmap name (including path)
0467 //  G4MAGNET::magfield_rescale = 1.;  // make consistent with expected Babar field strength of 1.4T
0468
0469   //---------------
0470   // Pythia Decayer
0471   //---------------
0472   // list of decay types in
0473   // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
0474   // default is All:
0475   // G4P6DECAYER::decayType = EDecayType::kAll;
0476
0477   // Initialize the selected subsystems
0478   G4Init();
0479
0480   //---------------------
0481   // GEANT4 Detector description
0482   //---------------------
0483   if (!Input::READHITS)
0484   {
0485     G4Setup();
0486   }
0487
0488   //------------------
0489   // Detector Division
0490   //------------------
0491
0492   if ((Enable::MBD && Enable::MBDRECO) || Enable::MBDFAKE) Mbd_Reco();
0493
0494   if (Enable::MVTX_CELL) Mvtx_Cells();
0495   if (Enable::INTT_CELL) Intt_Cells();
0496   if (Enable::TPC_CELL) TPC_Cells();
0497   if (Enable::MICROMEGAS_CELL) Micromegas_Cells();
0498
0499   if (Enable::CEMC_CELL) CEMC_Cells();
0500
0501   if (Enable::HCALIN_CELL) HCALInner_Cells();
0502
0503   if (Enable::HCALOUT_CELL) HCALOuter_Cells();
0504
0505   //-----------------------------
0506   // CEMC towering and clustering
0507   //-----------------------------
0508
0509   if (Enable::CEMC_TOWER) CEMC_Towers();
0510   if (Enable::CEMC_CLUSTER) CEMC_Clusters();
0511
0512   //--------------
0513   // EPD tile reconstruction
0514   //--------------
0515
0516   if (Enable::EPD_TILE) EPD_Tiles();
0517
0518   //-----------------------------
0519   // HCAL towering and clustering
0520   //-----------------------------
0521
0522   if (Enable::HCALIN_TOWER) HCALInner_Towers();
0523   if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
0524
0525   if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
0526   if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
0527
0528   // if enabled, do topoClustering early, upstream of any possible jet reconstruction
0529   if (Enable::TOPOCLUSTER) TopoClusterReco();
0530
0531   //--------------
0532   // SVTX tracking
0533   //--------------
0534   if(Enable::TRACKING_TRACK)
0535     {
0536       TrackingInit();
0537     }
0538   if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
0539   if (Enable::INTT_CLUSTER) Intt_Clustering();
0540   if (Enable::TPC_CLUSTER)
0541     {
0542       if(G4TPC::ENABLE_DIRECT_LASER_HITS || G4TPC::ENABLE_CENTRAL_MEMBRANE_HITS)
0543     {
0544       TPC_LaserClustering();
0545     }
0546       else
0547     {
0548       TPC_Clustering();
0549     }
0550     }
0551   if (Enable::MICROMEGAS_CLUSTER) Micromegas_Clustering();
0552
0553   if (Enable::TRACKING_TRACK)
0554   {
0555     Tracking_Reco();
0556   }
0557
0558
0559
0560   if(Enable::TRACKING_DIAGNOSTICS)
0561     {
0562       const std::string kshortFile = "./kshort_" + outputFile;
0563       const std::string residualsFile = "./residuals_" + outputFile;
0564
0565       G4KshortReconstruction(kshortFile);
0566       seedResiduals(residualsFile);
0567     }
0568
0569
0570   //-----------------
0571   // Global Vertexing
0572   //-----------------
0573
0574   if (Enable::GLOBAL_RECO && Enable::GLOBAL_FASTSIM)
0575   {
0576     cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
0577     gSystem->Exit(1);
0578   }
0579   if (Enable::GLOBAL_RECO)
0580   {
0581     Global_Reco();
0582   }
0583   else if (Enable::GLOBAL_FASTSIM)
0584   {
0585     Global_FastSim();
0586   }
0587
0588   //-----------------
0589   // Centrality Determination
0590   //-----------------
0591
0592   if (Enable::CENTRALITY)
0593   {
0594       Centrality();
0595   }
0596
0597   //-----------------
0598   // Calo Trigger Simulation
0599   //-----------------
0600
0601   if (Enable::CALOTRIGGER)
0602   {
0603     CaloTrigger_Sim();
0604   }
0605
0606   //---------
0607   // Jet reco
0608   //---------
0609
0610   if (Enable::JETS) Jet_Reco();
0611   if (Enable::HIJETS) HIJetReco();
0612
0613   if (Enable::PARTICLEFLOW) ParticleFlow();
0614
0615   //----------------------
0616   // Simulation evaluation
0617   //----------------------
0618   string outputroot = outputFile;
0619   string remove_this = ".root";
0620   size_t pos = outputroot.find(remove_this);
0621   if (pos != string::npos)
0622   {
0623     outputroot.erase(pos, remove_this.length());
0624   }
0625
0626   if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4svtx_eval.root");
0627
0628   if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
0629
0630   if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
0631
0632   if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
0633
0634   if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
0635
0636   if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");
0637
0638
0639
0640   if (Enable::USER) UserAnalysisInit();
0641  AnaTutorial *anaTutorial = new AnaTutorial("anaTutorial", outputroot + "_anaTutorial.root");
0642   anaTutorial->setMinJetPt(10.);
0643   anaTutorial->Verbosity(0);
0644   anaTutorial->analyzeTracks(true);
0645   anaTutorial->analyzeClusters(true);
0646   anaTutorial->analyzeJets(true);
0647   anaTutorial->analyzeTruth(false);
0648   se->registerSubsystem(anaTutorial);
0649   // Writes electrons from conversions to a new track map on the node tree
0650   // the ntuple file is for diagnostics, it is produced only if the flag is set in G4_Tracking.C
0651   if(G4TRACKING::filter_conversion_electrons) Filter_Conversion_Electrons(outputroot + "_secvert_ntuple.root");
0652
0653
0654   //======================
0655   // Run KFParticle on evt
0656   //======================
0657   if (Enable::KFPARTICLE && Input::UPSILON) KFParticle_Upsilon_Reco();
0658   if (Enable::KFPARTICLE && Input::DZERO) KFParticle_D0_Reco();
0659
0660   //----------------------
0661   // Standard QAs
0662   //----------------------
0663
0664   if (Enable::CEMC_QA) CEMC_QA();
0665   if (Enable::HCALIN_QA) HCALInner_QA();
0666   if (Enable::HCALOUT_QA) HCALOuter_QA();
0667
0668   if (Enable::JETS_QA) Jet_QA();
0669
0670   if (Enable::MVTX_QA) Mvtx_QA();
0671   if (Enable::INTT_QA) Intt_QA();
0672   if (Enable::TPC_QA) TPC_QA();
0673   if (Enable::MICROMEGAS_QA) Micromegas_QA();
0674   if (Enable::TRACKING_QA) Tracking_QA();
0675
0676   if (Enable::TRACKING_QA && Enable::CEMC_QA && Enable::HCALIN_QA && Enable::HCALOUT_QA) QA_G4CaloTracking();
0677
0678   if (Enable::TRACK_MATCHING) Track_Matching(outputroot + "_g4trackmatching.root");
0679
0680   //--------------
0681   // Set up Input Managers
0682   //--------------
0683
0684   InputManagers();
0685
0686   if (Enable::PRODUCTION)
0687   {
0688     Production_CreateOutputDir();
0689   }
0690
0691   if (Enable::DSTOUT)
0692   {
0693     string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
0694     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0695     if (Enable::DSTOUT_COMPRESS)
0696     {
0697       ShowerCompress();
0698       DstCompress(out);
0699     }
0700     se->registerOutputManager(out);
0701   }
0702   //-----------------
0703   // Event processing
0704   //-----------------
0705   if (Enable::DISPLAY)
0706   {
0707     DisplayOn();
0708
0709     gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
0710     gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
0711
0712     cout << "-------------------------------------------------" << endl;
0713     cout << "You are in event display mode. Run one event with" << endl;
0714     cout << "se->run(1)" << endl;
0715     cout << "Run Geant4 command with following examples" << endl;
0716     gROOT->ProcessLine("displaycmd()");
0717
0718     return 0;
0719   }
0720
0721   // if we use a negative number of events we go back to the command line here
0722   if (nEvents < 0)
0723   {
0724     return 0;
0725   }
0726   // if we run the particle generator and use 0 it'll run forever
0727   // for embedding it runs forever if the repeat flag is set
0728   if (nEvents == 0 && !Input::HEPMC && !Input::READHITS && INPUTEMBED::REPEAT)
0729   {
0730     cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0731     cout << "it will run forever, so I just return without running anything" << endl;
0732     return 0;
0733   }
0734
0735   se->skip(skip);
0736   se->run(nEvents);
0737   //  se->PrintTimer();
0738
0739   //-----
0740   // QA output
0741   //-----
0742
0743   if (Enable::QA) QA_Output(outputroot + "_qa.root");
0744
0745   //-----
0746   // Exit
0747   //-----
0748
0749 //  CDBInterface::instance()->Print(); // print used DB files
0750   se->End();
0751   std::cout << "All done" << std::endl;
0752   delete se;
0753   if (Enable::PRODUCTION)
0754   {
0755     Production_MoveOutput();
0756   }
0757
0758   gSystem->Exit(0);
0759   return 0;
0760 }
0761 #endif