AnaSvtxTracksForGenFit/macros/Fun4All_G4_sPHENIX_AnaGenFit.C

0001
0002 int Fun4All_G4_sPHENIX_AnaGenFit(
0003   const int nEvents = 10000,
0004   const char * inputFile = "/gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.5/single_particle/spacal1d/fieldmap/G4Hits_sPHENIX_e-_eta0_16GeV.root",
0005   const char * outputFile = "AnaSvtxTracksForGenFit.root"
0006 )
0007 {
0008   //===============
0009   // Input options
0010   //===============
0011
0012   // Either:
0013   // read previously generated g4-hits files, in this case it opens a DST and skips
0014   // the simulations step completely. The G4Setup macro is only loaded to get information
0015   // about the number of layers used for the cell reco code
0016   const bool readhits = false;
0017   // Or:
0018   // read files in HepMC format (typically output from event generators like hijing or pythia)
0019   const bool readhepmc = false; // read HepMC files
0020   // Or:
0021   // Use particle generator
0022   const bool runpythia8 = false;
0023   const bool runpythia6 = false;
0024
0025   //======================
0026   // What to run
0027   //======================
0028
0029   bool do_bbc = true;
0030
0031   bool do_pipe = true;
0032
0033   bool do_svtx = true;
0034   bool do_svtx_cell = true;
0035   bool do_svtx_track = true;
0036   bool do_svtx_eval = true;
0037
0038   bool do_preshower = false;
0039
0040   bool do_cemc = false;
0041   bool do_cemc_cell = false;
0042   bool do_cemc_twr = false;
0043   bool do_cemc_cluster = false;
0044   bool do_cemc_eval = false;
0045
0046   bool do_hcalin = false;
0047   bool do_hcalin_cell = false;
0048   bool do_hcalin_twr = false;
0049   bool do_hcalin_cluster = false;
0050   bool do_hcalin_eval = false;
0051
0052   bool do_magnet = false;
0053
0054   bool do_hcalout = false;
0055   bool do_hcalout_cell = false;
0056   bool do_hcalout_twr = false;
0057   bool do_hcalout_cluster = false;
0058   bool do_hcalout_eval = false;
0059
0060   bool do_global = false;
0061   bool do_global_fastsim = false;
0062
0063   bool do_jet_reco = false;
0064   bool do_jet_eval = false;
0065
0066   bool do_dst_compress = false;
0067
0068   //Option to convert DST to human command readable TTree for quick poke around the outputs
0069   bool do_DSTReader = false;
0070   //---------------
0071   // Load libraries
0072   //---------------
0073
0074   gSystem->Load("libfun4all.so");
0075   gSystem->Load("libg4detectors.so");
0076   gSystem->Load("libphhepmc.so");
0077   gSystem->Load("libg4testbench.so");
0078   gSystem->Load("libg4hough.so");
0079   gSystem->Load("libg4calo.so");
0080   gSystem->Load("libg4eval.so");
0081   gSystem->Load("libAnaSvtxTracksForGenFit.so");
0082
0083   // establish the geometry and reconstruction setup
0084   gROOT->LoadMacro("G4Setup_sPHENIX.C");
0085   G4Init(do_svtx, do_preshower, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe);
0086
0087   int absorberactive = 1; // set to 1 to make all absorbers active volumes
0088   //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0089   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)
0090   const float magfield_rescale = 1.4 / 1.5; // scale the map to a 1.4 T field
0091
0092   //---------------
0093   // Fun4All server
0094   //---------------
0095
0096   Fun4AllServer *se = Fun4AllServer::instance();
0097   se->Verbosity(0);
0098   // just if we set some flags somewhere in this macro
0099   recoConsts *rc = recoConsts::instance();
0100   // By default every random number generator uses
0101   // PHRandomSeed() which reads /dev/urandom to get its seed
0102   // if the RANDOMSEED flag is set its value is taken as seed
0103   // You ca neither set this to a random value using PHRandomSeed()
0104   // which will make all seeds identical (not sure what the point of
0105   // this would be:
0106   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0107   // or set it to a fixed value so you can debug your code
0108   // rc->set_IntFlag("RANDOMSEED", 12345);
0109
0110   //-----------------
0111   // Event generation
0112   //-----------------
0113
0114   if (readhits)
0115   {
0116     // Get the hits from a file
0117     // The input manager is declared later
0118   }
0119   else if (readhepmc)
0120   {
0121     // this module is needed to read the HepMC records into our G4 sims
0122     // but only if you read HepMC input files
0123     HepMCNodeReader *hr = new HepMCNodeReader();
0124     se->registerSubsystem(hr);
0125   }
0126   else if (runpythia8)
0127   {
0128     gSystem->Load("libPHPythia8.so");
0129
0130     PHPythia8* pythia8 = new PHPythia8();
0131     // see coresoftware/generators/PHPythia8 for example config
0132     pythia8->set_config_file("phpythia8.cfg");
0133     se->registerSubsystem(pythia8);
0134
0135     HepMCNodeReader *hr = new HepMCNodeReader();
0136     se->registerSubsystem(hr);
0137   }
0138   else if (runpythia6)
0139   {
0140     gSystem->Load("libPHPythia6.so");
0141
0142     PHPythia6 *pythia6 = new PHPythia6();
0143     pythia6->set_config_file("phpythia6.cfg");
0144     se->registerSubsystem(pythia6);
0145
0146     HepMCNodeReader *hr = new HepMCNodeReader();
0147     se->registerSubsystem(hr);
0148   }
0149   else
0150   {
0151     // toss low multiplicity dummy events
0152     PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0153     gen->add_particles("mu+", 1); // mu+,e+,proton,pi+,Upsilon
0154     // gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
0155     if (readhepmc) {
0156       gen->set_reuse_existing_vertex(true);
0157       gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0158     } else {
0159       gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0160                                             PHG4SimpleEventGenerator::Uniform,
0161                                             PHG4SimpleEventGenerator::Uniform);
0162       gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
0163       gen->set_vertex_distribution_width(0.0, 0.0, 5.0);
0164     }
0165     gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0166     gen->set_vertex_size_parameters(0.0, 0.0);
0167     gen->set_eta_range(-0.5, 0.5);
0168     gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0169     gen->set_pt_range(0.1, 10.0);
0170     gen->Embed(1);
0171     gen->Verbosity(0);
0172     se->registerSubsystem(gen);
0173   }
0174
0175   if (!readhits)
0176   {
0177     //---------------------
0178     // Detector description
0179     //---------------------
0180
0181     G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
0182             do_svtx, do_preshower, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe, magfield_rescale);
0183   }
0184
0185   //---------
0186   // BBC Reco
0187   //---------
0188
0189   if (do_bbc)
0190   {
0191     gROOT->LoadMacro("G4_Bbc.C");
0192     BbcInit();
0193     Bbc_Reco();
0194   }
0195   //------------------
0196   // Detector Division
0197   //------------------
0198
0199   if (do_svtx_cell) Svtx_Cells();
0200
0201   if (do_cemc_cell) CEMC_Cells();
0202
0203   if (do_hcalin_cell) HCALInner_Cells();
0204
0205   if (do_hcalout_cell) HCALOuter_Cells();
0206
0207   //-----------------------------
0208   // CEMC towering and clustering
0209   //-----------------------------
0210
0211   if (do_cemc_twr) CEMC_Towers();
0212   if (do_cemc_cluster) CEMC_Clusters();
0213
0214   //-----------------------------
0215   // HCAL towering and clustering
0216   //-----------------------------
0217
0218   if (do_hcalin_twr) HCALInner_Towers();
0219   if (do_hcalin_cluster) HCALInner_Clusters();
0220
0221   if (do_hcalout_twr) HCALOuter_Towers();
0222   if (do_hcalout_cluster) HCALOuter_Clusters();
0223
0224   if (do_dst_compress) ShowerCompress();
0225
0226   //--------------
0227   // SVTX tracking
0228   //--------------
0229
0230   if (do_svtx_track) Svtx_Reco();
0231
0232   //-----------------
0233   // Global Vertexing
0234   //-----------------
0235
0236   if (do_global)
0237   {
0238     gROOT->LoadMacro("G4_Global.C");
0239     Global_Reco();
0240   }
0241
0242   else if (do_global_fastsim)
0243   {
0244     gROOT->LoadMacro("G4_Global.C");
0245     Global_FastSim();
0246   }
0247
0248   //---------
0249   // Jet reco
0250   //---------
0251
0252   if (do_jet_reco)
0253   {
0254     gROOT->LoadMacro("G4_Jets.C");
0255     Jet_Reco();
0256   }
0257   //----------------------
0258   // Simulation evaluation
0259   //----------------------
0260
0261   if (do_svtx_eval) Svtx_Eval("g4svtx_eval.root");
0262
0263   if (do_cemc_eval) CEMC_Eval("g4cemc_eval.root");
0264
0265   if (do_hcalin_eval) HCALInner_Eval("g4hcalin_eval.root");
0266
0267   if (do_hcalout_eval) HCALOuter_Eval("g4hcalout_eval.root");
0268
0269   if (do_jet_eval) Jet_Eval("g4jet_eval.root");
0270
0271   //--------------
0272   // IO management
0273   //--------------
0274
0275   if (readhits)
0276   {
0277     // Hits file
0278     Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
0279     hitsin->fileopen(inputFile);
0280     se->registerInputManager(hitsin);
0281   }
0282   if (readhepmc)
0283   {
0284     Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0285     se->registerInputManager( in );
0286     se->fileopen( in->Name().c_str(), inputFile );
0287   }
0288   else
0289   {
0290     // for single particle generators we just need something which drives
0291     // the event loop, the Dummy Input Mgr does just that
0292     Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0293     se->registerInputManager( in );
0294   }
0295
0296   if (do_DSTReader)
0297   {
0298     //Convert DST to human command readable TTree for quick poke around the outputs
0299     gROOT->LoadMacro("G4_DSTReader.C");
0300
0301     G4DSTreader( outputFile, //
0302                  /*int*/ absorberactive ,
0303                  /*bool*/ do_svtx ,
0304                  /*bool*/ do_preshower ,
0305                  /*bool*/ do_cemc ,
0306                  /*bool*/ do_hcalin ,
0307                  /*bool*/ do_magnet ,
0308                  /*bool*/ do_hcalout ,
0309                  /*bool*/ do_cemc_twr ,
0310                  /*bool*/ do_hcalin_twr ,
0311                  /*bool*/ do_magnet  ,
0312                  /*bool*/ do_hcalout_twr
0313                );
0314   }
0315
0316   // Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0317   // if (do_dst_compress) DstCompress(out);
0318   // se->registerOutputManager(out);
0319
0320   //-----------------
0321   // Analysis Module
0322   //-----------------
0323
0324   AnaSvtxTracksForGenFit *ana = new AnaSvtxTracksForGenFit();
0325   ana->set_filename( outputFile );
0326   se->registerSubsystem( ana );
0327
0328   //-----------------
0329   // Event processing
0330   //-----------------
0331   if (nEvents < 0)
0332   {
0333     return;
0334   }
0335   // if we run the particle generator and use 0 it'll run forever
0336   if (nEvents == 0 && !readhits && !readhepmc)
0337   {
0338     cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0339     cout << "it will run forever, so I just return without running anything" << endl;
0340     return;
0341   }
0342
0343   se->run(nEvents);
0344
0345   //-----
0346   // Exit
0347   //-----
0348
0349   se->End();
0350   std::cout << "All done" << std::endl;
0351   delete se;
0352   gSystem->Exit(0);
0353 }