pythiaInfo/macro/Fun4AllPythia.C

0001 #include <fun4all/Fun4AllDstOutputManager.h>
0002 #include <fun4all/Fun4AllOutputManager.h>
0003 #include <fun4all/Fun4AllServer.h>
0004 #include <fun4all/Fun4AllSyncManager.h>
0005 #include <fun4all/Fun4AllUtils.h>
0006
0007
0008 #include <phool/PHRandomSeed.h>
0009 #include <phool/recoConsts.h>
0010
0011 #include <stdlib.h>
0012
0013 #include <g4main/PHG4Reco.h>
0014 #include <g4main/PHG4TruthSubsystem.h>
0015
0016 #include <GlobalVariables.C>
0017
0018 #include <G4_Input.C>
0019 #include <phpythia8/PHPy8JetTrigger.h>
0020
0021 #include <jethistogrammer/jetHistogrammer.h>
0022
0023 #include <g4jets/FastJetAlgo.h>
0024 #include <g4jets/JetReco.h>
0025 #include <g4jets/TruthJetInput.h>
0026 #include <time.h>
0027
0028 R__LOAD_LIBRARY(libjetHistogrammer.so)
0029
0030 void Fun4AllPythia(int nEvents = 1,
0031            const string &Jet_Trigger = "Jet10", // or "PhotonJet"
0032            int doCrossSection = 1,
0033            const string &outname = "out"
0034 )
0035 {
0036
0037   clock_t tStart = clock();
0038
0039
0040   Fun4AllServer *se = Fun4AllServer::instance();
0041   se -> Verbosity(0);
0042
0043   Input::PYTHIA8 = true;
0044   Input::VERBOSITY = 0;
0045
0046
0047   if (Jet_Trigger == "PhotonJet")
0048     {
0049       PYTHIA8::config_file = "phpythia8_JS_GJ_MDC2.cfg";
0050     }
0051   else if (Jet_Trigger == "Jet10")
0052     {
0053       PYTHIA8::config_file = "phpythia8_15GeV_JS_MDC2.cfg";
0054     }
0055   else if (Jet_Trigger == "Jet30")
0056     {
0057       PYTHIA8::config_file = "phpythia8_30GeV_JS_MDC2.cfg";
0058     }
0059   else if(Jet_Trigger == "MB")
0060     {
0061       PYTHIA8::config_file = "phpythia8_minBias_MDC2.cfg";
0062     }
0063   else
0064     {
0065       std::cout << "Invalid jet trigger " << Jet_Trigger << std::endl;
0066       gSystem->Exit(1);
0067     }
0068
0069   InputInit();
0070
0071   if(Jet_Trigger != "MB")
0072     {
0073       PHPy8JetTrigger *p8_js_signal_trigger = new PHPy8JetTrigger();
0074       p8_js_signal_trigger->SetEtaHighLow(1.5,-1.5); // Set eta acceptance for particles into the jet between +/- 1.5
0075       p8_js_signal_trigger->SetJetR(0.4);      //Set the radius for the trigger jet
0076       p8_js_signal_trigger->Verbosity(0);;
0077       if (Jet_Trigger == "Jet10")
0078     {
0079       p8_js_signal_trigger->SetMinJetPt(10); // require a 10 GeV minimum pT jet in the event
0080     }
0081       else if (Jet_Trigger == "Jet30")
0082     {
0083       p8_js_signal_trigger->SetMinJetPt(30); // require a 30 GeV minimum pT jet in the event
0084     }
0085       else if (Jet_Trigger == "PhotonJet")
0086     {
0087       delete p8_js_signal_trigger;
0088       p8_js_signal_trigger = nullptr;
0089       cout << "no cut for PhotonJet" << endl;
0090     }
0091       else
0092     {
0093       cout << "invalid jet trigger: " << Jet_Trigger << endl;
0094       gSystem->Exit(1);
0095     }
0096
0097       if (p8_js_signal_trigger)
0098     {
0099       INPUTGENERATOR::Pythia8->register_trigger(p8_js_signal_trigger);
0100       INPUTGENERATOR::Pythia8->set_trigger_AND();
0101     }
0102     }
0103   Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
0104   INPUTGENERATOR::Pythia8 -> Verbosity(doCrossSection*10);
0105
0106   InputRegister();
0107   InputManagers();
0108   PHG4Reco *reco = new PHG4Reco();
0109   reco -> set_field(0);
0110   reco -> SetWorldMaterial("G4_Galactic");
0111   reco -> SetWorldSizeX(100);
0112   reco -> SetWorldSizeY(100);
0113   reco -> SetWorldSizeZ(100);
0114   reco -> save_DST_geometry(false);
0115
0116   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0117   reco -> registerSubsystem(truth);
0118   se -> registerSubsystem(reco);
0119
0120   JetReco *truthJetReco = new JetReco();
0121   TruthJetInput *jetInput = new TruthJetInput(Jet::PARTICLE);
0122   truthJetReco -> add_input(jetInput);
0123   truthJetReco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.4), "AntiKt_Truth_r04");
0124   truthJetReco->set_algo_node("ANTIKT");
0125   truthJetReco->set_input_node("TRUTH");
0126   truthJetReco->Verbosity(0);
0127
0128   se -> registerSubsystem(truthJetReco);
0129
0130   string out = outname + Jet_Trigger + ".root";
0131   jetHistogrammer *jetHistos = new jetHistogrammer("jetHistogrammer",out.c_str());
0132   se -> registerSubsystem(jetHistos);
0133
0134
0135   se -> run(nEvents);
0136
0137   se -> End();
0138   std::cout << "All done" << std::endl;
0139
0140
0141   delete se;
0142   std::cout << "Total runtime: " << double(clock() - tStart) / (double)CLOCKS_PER_SEC << std::endl;;
0143
0144   gSystem -> Exit(0);
0145   return 0;
0146 }
0147
0148