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