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