PhotonJet/example_macros/Fun4All_G4_fsPHENIX_pi0.C

0001
0002 int Fun4All_G4_fsPHENIX_pi0(
0003                const int nEvents = 10,
0004            const char * inputFile = "/sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/fieldmap/G4Hits_sPHENIX_e-_eta0_8GeV-0002.root",
0005                const char * outputFile = "G4fsPHENIX.root",
0006            const char * embed_input_file = "/sphenix/sim/sim01/production/2016-07-12/sHijing/spacal2d/G4Hits_sPHENIX_sHijing-0-4.4fm.list"
0007                )
0008 {
0009   // Set the number of TPC layer
0010   const int n_TPC_layers = 40;  // use 60 for backward compatibility only
0011
0012   //===============
0013   // Input options
0014   //===============
0015
0016   // Either:
0017   // read previously generated g4-hits files, in this case it opens a DST and skips
0018   // the simulations step completely. The G4Setup macro is only loaded to get information
0019   // about the number of layers used for the cell reco code
0020   const bool readhits = false;
0021   // Or:
0022   // read files in HepMC format (typically output from event generators like hijing or pythia)
0023   const bool readhepmc = false; // read HepMC files
0024   // Or:
0025   // Use particle generator
0026   const bool runpythia8 = false;
0027   const bool runpythia6 = false;
0028   // And
0029   // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
0030   // 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
0031   // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0032   const bool do_embedding = false;
0033
0034   //======================
0035   // What to run
0036   //======================
0037
0038   bool do_bbc = true;
0039
0040   bool do_pipe = true;
0041
0042   bool do_svtx = true;
0043   bool do_svtx_cell = do_svtx && true;
0044   bool do_svtx_track = do_svtx_cell && true;
0045   bool do_svtx_eval = do_svtx_track && false;
0046
0047   bool do_cemc = true;
0048   bool do_cemc_cell = do_cemc && true;
0049   bool do_cemc_twr = do_cemc_cell && true;
0050   bool do_cemc_cluster = do_cemc_twr && true;
0051   bool do_cemc_eval = do_cemc_cluster && false;
0052
0053   bool do_hcalin = true;
0054   bool do_hcalin_cell = do_hcalin && true;
0055   bool do_hcalin_twr = do_hcalin_cell && true;
0056   bool do_hcalin_cluster = do_hcalin_twr && true;
0057   bool do_hcalin_eval = do_hcalin_cluster && false;
0058
0059   bool do_magnet = true;
0060
0061   bool do_hcalout = true;
0062   bool do_hcalout_cell = do_hcalout && true;
0063   bool do_hcalout_twr = do_hcalout_cell && true;
0064   bool do_hcalout_cluster = do_hcalout_twr && true;
0065   bool do_hcalout_eval = do_hcalout_cluster && false;
0066
0067   bool do_global = true;
0068   bool do_global_fastsim = false;
0069
0070   bool do_jet_reco = false;
0071   bool do_jet_eval = do_jet_reco && true;
0072
0073   bool do_fwd_jet_reco = false;
0074   bool do_fwd_jet_eval = do_fwd_jet_reco && false;
0075
0076   // fsPHENIX geometry
0077
0078   bool do_FGEM = true;
0079   bool do_FGEM_track = do_FGEM &&  false;
0080
0081   bool do_FEMC = true;
0082   bool do_FEMC_cell = do_FEMC && true;
0083   bool do_FEMC_twr = do_FEMC_cell && true;
0084   bool do_FEMC_cluster = do_FEMC_twr && true;
0085
0086   bool do_FHCAL = true;
0087   bool do_FHCAL_cell = do_FHCAL && true;
0088   bool do_FHCAL_twr = do_FHCAL_cell && true;
0089   bool do_FHCAL_cluster = do_FHCAL_twr && true;
0090
0091   bool do_dst_compress = false;
0092
0093   //Option to convert DST to human command readable TTree for quick poke around the outputs
0094   bool do_DSTReader = true;
0095   //---------------
0096   // Load libraries
0097   //---------------
0098
0099   gSystem->Load("libfun4all.so");
0100   gSystem->Load("libg4detectors.so");
0101   gSystem->Load("libphhepmc.so");
0102   gSystem->Load("libg4testbench.so");
0103   gSystem->Load("libg4hough.so");
0104   gSystem->Load("libg4eval.so");
0105
0106   // establish the geometry and reconstruction setup
0107   gROOT->LoadMacro("G4Setup_fsPHENIX.C");
0108   G4Init(do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,do_FGEM,do_FEMC,do_FHCAL,n_TPC_layers);
0109
0110   int absorberactive = 0; // set to 1 to make all absorbers active volumes
0111   //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0112   const string magfield = "/phenix/upgrades/decadal/fieldmaps/fsPHENIX.2d.root"; // fsPHENIX field map by Cesar Luiz da Silva <slash@rcf.rhic.bnl.gov>, sPHENIX + piston
0113   const float magfield_rescale = 1.0; // already adjusted to 1.4T central field
0114
0115   //---------------
0116   // Fun4All server
0117   //---------------
0118
0119   Fun4AllServer *se = Fun4AllServer::instance();
0120 //  se->Verbosity(0); // uncomment for batch production running with minimal output messages
0121   se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running
0122   // just if we set some flags somewhere in this macro
0123   recoConsts *rc = recoConsts::instance();
0124   // By default every random number generator uses
0125   // PHRandomSeed() which reads /dev/urandom to get its seed
0126   // if the RANDOMSEED flag is set its value is taken as seed
0127   // You can either set this to a random value using PHRandomSeed()
0128   // which will make all seeds identical (not sure what the point of
0129   // this would be:
0130   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0131   // or set it to a fixed value so you can debug your code
0132   // rc->set_IntFlag("RANDOMSEED", 12345);
0133
0134   //-----------------
0135   // Event generation
0136   //-----------------
0137
0138   if (readhits)
0139     {
0140       // Get the hits from a file
0141       // The input manager is declared later
0142       if (do_embedding)
0143        {
0144          cout <<"Do not support read hits and embed background at the same time."<<endl;
0145          exit(1);
0146        }
0147     }
0148   else if (readhepmc)
0149     {
0150     // action is performed in later stage at the input manager level
0151     }
0152   else if (runpythia8)
0153     {
0154       gSystem->Load("libPHPythia8.so");
0155
0156       PHPythia8* pythia8 = new PHPythia8();
0157       // see coresoftware/generators/PHPythia8 for example config
0158       pythia8->set_config_file("phpythia8.cfg");
0159       se->registerSubsystem(pythia8);
0160     }
0161   else if (runpythia6)
0162     {
0163       gSystem->Load("libPHPythia6.so");
0164
0165       PHPythia6 *pythia6 = new PHPythia6();
0166       pythia6->set_config_file("phpythia6.cfg");
0167       se->registerSubsystem(pythia6);
0168     }
0169   else
0170     {
0171       // toss low multiplicity dummy events
0172       PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0173       //gen->add_particles("e-",5); // mu+,e+,proton,pi+,Upsilon
0174       //gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
0175       //gen->add_particles("pi-",1); // mu-,e-,anti_proton,pi-
0176       gen->add_particles(111,1);
0177
0178       if (readhepmc || do_embedding) {
0179     gen->set_reuse_existing_vertex(true);
0180     gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
0181       } else {
0182     gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0183                            PHG4SimpleEventGenerator::Uniform,
0184                            PHG4SimpleEventGenerator::Uniform);
0185     gen->set_vertex_distribution_mean(0.0,0.0,0.0);
0186     gen->set_vertex_distribution_width(0.0,0.0,5.0);
0187       }
0188       gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0189       gen->set_vertex_size_parameters(0.0,0.0);
0190       gen->set_eta_range(2.2, 2.5);
0191       //gen->set_eta_range(3.0, 3.0); //fsPHENIX FWD
0192       gen->set_phi_range(-1.0*TMath::Pi(), 1.0*TMath::Pi());
0193       //gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
0194       gen->set_p_range(5.0, 8.0);
0195       //gen->Embed(1);
0196       gen->Verbosity(0);
0197       se->registerSubsystem(gen);
0198     }
0199
0200   if (!readhits)
0201     {
0202       //---------------------
0203       // Detector description
0204       //---------------------
0205
0206       G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
0207           do_svtx, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe,
0208           do_FGEM, do_FEMC, do_FHCAL,
0209           magfield_rescale);
0210
0211     }
0212
0213   //---------
0214   // BBC Reco
0215   //---------
0216
0217   if (do_bbc)
0218     {
0219       gROOT->LoadMacro("G4_Bbc.C");
0220       BbcInit();
0221       Bbc_Reco();
0222     }
0223
0224   //------------------
0225   // Detector Division
0226   //------------------
0227
0228   if (do_svtx_cell) Svtx_Cells();
0229
0230   if (do_cemc_cell) CEMC_Cells();
0231
0232   if (do_hcalin_cell) HCALInner_Cells();
0233
0234   if (do_hcalout_cell) HCALOuter_Cells();
0235
0236   if (do_FEMC_cell) FEMC_Cells();
0237   if (do_FHCAL_cell) FHCAL_Cells();
0238
0239   //-----------------------------
0240   // CEMC towering and clustering
0241   //-----------------------------
0242
0243   if (do_cemc_twr) CEMC_Towers();
0244   if (do_cemc_cluster) CEMC_Clusters();
0245
0246   //-----------------------------
0247   // HCAL towering and clustering
0248   //-----------------------------
0249
0250   if (do_hcalin_twr) HCALInner_Towers();
0251   if (do_hcalin_cluster) HCALInner_Clusters();
0252
0253   if (do_hcalout_twr) HCALOuter_Towers();
0254   if (do_hcalout_cluster) HCALOuter_Clusters();
0255
0256   if (do_FEMC_twr) FEMC_Towers();
0257   if (do_FEMC_cluster) FEMC_Clusters();
0258
0259   if (do_FHCAL_twr) FHCAL_Towers();
0260   if (do_FHCAL_cluster) FHCAL_Clusters();
0261
0262   if (do_dst_compress) ShowerCompress();
0263
0264   //--------------
0265   // SVTX tracking
0266   //--------------
0267
0268   if (do_svtx_track) Svtx_Reco();
0269
0270   //--------------
0271   // FGEM tracking
0272   //--------------
0273
0274   if(do_FGEM_track) FGEM_FastSim_Reco();
0275
0276   //-----------------
0277   // Global Vertexing
0278   //-----------------
0279
0280   if (do_global)
0281     {
0282       gROOT->LoadMacro("G4_Global.C");
0283       Global_Reco();
0284     }
0285
0286   else if (do_global_fastsim)
0287     {
0288       gROOT->LoadMacro("G4_Global.C");
0289       Global_FastSim();
0290     }
0291
0292   //---------
0293   // Jet reco
0294   //---------
0295
0296   if (do_jet_reco)
0297     {
0298       gROOT->LoadMacro("G4_Jets.C");
0299       Jet_Reco();
0300     }
0301
0302   if (do_fwd_jet_reco)
0303     {
0304       gROOT->LoadMacro("G4_FwdJets.C");
0305       Jet_FwdReco();
0306     }
0307   //----------------------
0308   // Simulation evaluation
0309   //----------------------
0310
0311   if (do_svtx_eval) Svtx_Eval("g4svtx_eval.root");
0312
0313   if (do_cemc_eval) CEMC_Eval("g4cemc_eval.root");
0314
0315   if (do_hcalin_eval) HCALInner_Eval("g4hcalin_eval.root");
0316
0317   if (do_hcalout_eval) HCALOuter_Eval("g4hcalout_eval.root");
0318
0319   if (do_jet_eval) Jet_Eval("g4jet_eval.root");
0320
0321   if (do_fwd_jet_eval) Jet_FwdEval("g4fwdjet_eval.root");
0322
0323
0324   gSystem->Load("libPhotons.so");
0325   Forward_pi0s *photons = new Forward_pi0s(outputFile);
0326   photons->Verbosity(0);
0327
0328   photons->set_cluspt_mincut(0.5);
0329   photons->set_eta_lowhigh(1.1,3.5);
0330   se->registerSubsystem(photons);
0331
0332
0333
0334
0335   //--------------
0336   // IO management
0337   //--------------
0338
0339   if (readhits)
0340     {
0341       // Hits file
0342       Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
0343       hitsin->fileopen(inputFile);
0344       se->registerInputManager(hitsin);
0345     }
0346   if (do_embedding)
0347     {
0348       if (embed_input_file == NULL)
0349         {
0350           cout << "Missing embed_input_file! Exit";
0351           exit(3);
0352         }
0353
0354       Fun4AllDstInputManager *in1 = new Fun4AllNoSyncDstInputManager("DSTinEmbed");
0355       //      in1->AddFile(embed_input_file); // if one use a single input file
0356       in1->AddListFile(embed_input_file); // RecommendedL: if one use a text list of many input files
0357       se->registerInputManager(in1);
0358     }
0359   if (readhepmc)
0360     {
0361       Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0362       se->registerInputManager( in );
0363       se->fileopen( in->Name().c_str(), inputFile );
0364     }
0365   else
0366     {
0367       // for single particle generators we just need something which drives
0368       // the event loop, the Dummy Input Mgr does just that
0369       Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0370       se->registerInputManager( in );
0371     }
0372
0373   if (do_DSTReader)
0374     {
0375       //Convert DST to human command readable TTree for quick poke around the outputs
0376       gROOT->LoadMacro("G4_DSTReader_fsPHENIX.C");
0377
0378       G4DSTreader_fsPHENIX( outputFile, //
0379           /*int*/ absorberactive ,
0380           /*bool*/ do_svtx ,
0381           /*bool*/ do_cemc ,
0382           /*bool*/ do_hcalin ,
0383           /*bool*/ do_magnet ,
0384           /*bool*/ do_hcalout ,
0385           /*bool*/ do_cemc_twr ,
0386           /*bool*/ do_hcalin_twr ,
0387           /*bool*/ do_magnet  ,
0388       /*bool*/ do_hcalout_twr,
0389     /*bool*/ do_FGEM,
0390       /*bool*/ do_FHCAL,
0391       /*bool*/ do_FHCAL_twr,
0392       /*bool*/ do_FEMC,
0393       /*bool*/ do_FEMC_twr
0394           );
0395     }
0396
0397   //Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0398   //if (do_dst_compress) DstCompress(out);
0399   //se->registerOutputManager(out);
0400
0401   if (nEvents == 0 && !readhits && !readhepmc)
0402     {
0403       cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0404       cout << "it will run forever, so I just return without running anything" << endl;
0405       return;
0406     }
0407
0408   if (nEvents < 0)
0409     {
0410       PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0411       g4->ApplyCommand("/control/execute vis.mac");
0412       //g4->StartGui();
0413       se->run(1);
0414
0415       se->End();
0416       std::cout << "All done" << std::endl;
0417
0418
0419       std::cout << "==== Useful display commands ==" << std::endl;
0420       cout << "draw axis: " << endl;
0421       cout << " G4Cmd(\"/vis/scene/add/axes 0 0 0 50 cm\")" << endl;
0422       cout << "zoom" << endl;
0423       cout << " G4Cmd(\"/vis/viewer/zoom 1\")" << endl;
0424       cout << "viewpoint:" << endl;
0425       cout << " G4Cmd(\"/vis/viewer/set/viewpointThetaPhi 0 0\")" << endl;
0426       cout << "panTo:" << endl;
0427       cout << " G4Cmd(\"/vis/viewer/panTo 0 0 cm\")" << endl;
0428       cout << "print to eps:" << endl;
0429       cout << " G4Cmd(\"/vis/ogl/printEPS\")" << endl;
0430       cout << "set background color:" << endl;
0431       cout << " G4Cmd(\"/vis/viewer/set/background white\")" << endl;
0432       std::cout << "===============================" << std::endl;
0433     }
0434   else
0435     {
0436
0437       se->run(nEvents);
0438
0439       se->End();
0440       std::cout << "All done" << std::endl;
0441       delete se;
0442       gSystem->Exit(0);
0443     }
0444
0445 }
0446
0447
0448 void
0449 G4Cmd(const char * cmd)
0450 {
0451   Fun4AllServer *se = Fun4AllServer::instance();
0452   PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0453   g4->ApplyCommand(cmd);
0454 }