CalorimeterJets/macro/G4_HIJetReco.C

0001 #ifndef MACRO_G4HIJETRECO_C
0002 #define MACRO_G4HIJETRECO_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <g4jets/FastJetAlgo.h>
0007 #include <g4jets/JetReco.h>
0008 #include <g4jets/TowerJetInput.h>
0009 #include <g4jets/TruthJetInput.h>
0010 #include <g4centrality/PHG4CentralityReco.h>
0011
0012 #include <jetbackground/CopyAndSubtractJets.h>
0013 #include <jetbackground/DetermineTowerBackground.h>
0014 #include <jetbackground/FastJetAlgoSub.h>
0015 #include <jetbackground/RetowerCEMC.h>
0016 #include <jetbackground/SubtractTowers.h>
0017 #include <jetbackground/SubtractTowersCS.h>
0018
0019 #include <fun4all/Fun4AllServer.h>
0020 #include <fun4all/SubsysReco.h>
0021 #include <fun4all/Fun4AllInputManager.h>
0022 #include <fun4all/Fun4AllDstInputManager.h>
0023
0024 #include <fun4all/Fun4AllDstOutputManager.h>
0025 #include <fun4all/Fun4AllOutputManager.h>
0026
0027 #include <phool/PHRandomSeed.h>
0028 #include <phool/recoConsts.h>
0029
0030 #include <myjetanalysis/MyJetAnalysis.h>
0031 #include <negativetowerremover/NegativeTowerRemover.h>
0032
0033 #include <fastjet/PseudoJet.hh>
0034
0035 R__LOAD_LIBRARY(libfun4all.so)
0036 R__LOAD_LIBRARY(libg4jets.so)
0037 R__LOAD_LIBRARY(libjetbackground.so)
0038 R__LOAD_LIBRARY(libmyjetanalysis.so)
0039 R__LOAD_LIBRARY(libg4centrality.so)
0040
0041 namespace Enable
0042 {
0043   bool HIJETS = false;
0044   int HIJETS_VERBOSITY = 0;
0045 }  // namespace Enable
0046
0047 namespace G4HIJETS
0048 {
0049   bool do_flow = false;
0050   bool do_CS = false;
0051 }  // namespace G4HIJETS
0052
0053 void G4_HIJetReco(const char *filelisttruth = "dst_truth.list",
0054                   const char *filelisttruthjet = "dst_truth_jet.list",
0055                   const char *filelisthits = "dst_bbc_g4hit.list",
0056           const char *filelistcalo = "dst_calo_cluster.list",
0057           const string &outname = "output_000002",
0058                   int n_skip = 2,
0059                   int n_event = 5
0060              )
0061    {
0062    gSystem->Load("libg4dst");
0063    gSystem->Load("libmyjetanalysis");
0064
0065   int verbosity = 0;//std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY);
0066
0067   //---------------
0068   // Fun4All server
0069   //--------------1-
0070
0071   Fun4AllServer *se = Fun4AllServer::instance();
0072   se->Verbosity(verbosity);
0073
0074   PHG4CentralityReco *cent = new PHG4CentralityReco();
0075   cent->Verbosity(0);
0076   cent->GetCalibrationParameters().ReadFromFile("centrality", "xml", 0, 0, string(getenv("CALIBRATIONROOT")) + string("/Centrality/"));
0077   se->registerSubsystem( cent );
0078
0079 //JetRecolysis::InitRun" << end`l;
0080 //  //m_jetEvalStack = shared_ptr<JetEvalStack>(new JetEvalStack(topNode, m_recoJetName, m_truthJetName));
0081 //    //m_jetEvalStack->get_stvx_eval_stack()->set_use_initial_vertex(initial_vertex);
0082 //      return Fun4AllReturnCodes::EVENT_OK;
0083 //      *truthjetreco = new JetReco();
0084   //TruthJetInput *tji = new TruthJetInput(Jet::PARTICLE);
0085   //tji->add_embedding_flag(1);  // changes depending on signal vs. embedded
0086   //truthjetreco->add_input(tji);
0087   //truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.2), "AntiKt_Truth_r02_sub");
0088   //truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.3), "AntiKt_Truth_r03_sub");
0089   //truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.4), "AntiKt_Truth_r04_sub");
0090   //truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.5), "AntiKt_Truth_r05_sub");
0091   ///truthjetreco->set_algo_node("ANTIKT");
0092   //truthjetreco->set_input_node("TRUTH");
0093   //truthjetreco->Verbosity(verbosity);
0094   //se->registerSubsystem(truthjetreco);
0095
0096   RetowerCEMC *rcemc = new RetowerCEMC();
0097   rcemc->Verbosity(verbosity);
0098   se->registerSubsystem(rcemc);
0099
0100   JetReco *towerjetreco = new JetReco();
0101   //towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWER));
0102   towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWER_RETOWER));
0103   towerjetreco->add_input(new TowerJetInput(Jet::HCALIN_TOWER));
0104   towerjetreco->add_input(new TowerJetInput(Jet::HCALOUT_TOWER));
0105   towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.2), "AntiKt_Tower_HIRecoSeedsRaw_r02");
0106   towerjetreco->set_algo_node("ANTIKT");
0107   towerjetreco->set_input_node("TOWER");
0108   towerjetreco->Verbosity(verbosity);
0109   se->registerSubsystem(towerjetreco);
0110
0111   DetermineTowerBackground *dtb = new DetermineTowerBackground();
0112   dtb->SetBackgroundOutputName("TowerBackground_Sub1");
0113   dtb->SetFlow(G4HIJETS::do_flow);
0114   dtb->SetSeedType(0);
0115   dtb->SetSeedJetD(3);
0116   dtb->Verbosity(verbosity);
0117   se->registerSubsystem(dtb);
0118
0119   CopyAndSubtractJets *casj = new CopyAndSubtractJets();
0120   casj->SetFlowModulation(G4HIJETS::do_flow);
0121   //casj->SetBackgroundInputName("TowerBackground_Sub1");
0122   //casj->SetJetInputName("AntiKt_Tower_HIRecoSeedsRaw_r02");
0123   //casj->SetJetOutputName("AntiKt_Tower_HIRecoSeedsSub_r02");
0124   casj->Verbosity(verbosity);
0125   se->registerSubsystem(casj);
0126
0127   DetermineTowerBackground *dtb2 = new DetermineTowerBackground();
0128   dtb2->SetBackgroundOutputName("TowerBackground_Sub2");
0129   dtb2->SetFlow(G4HIJETS::do_flow);
0130   dtb2->SetSeedType(1);
0131   dtb2->SetSeedJetPt(7);
0132   dtb2->Verbosity(verbosity);
0133   se->registerSubsystem(dtb2);
0134
0135   //CopyAndSubtractJets *casj_final = new CopyAndSubtractJets();
0136   //casj_final->SetFlowModulation(G4HIJETS::do_flow);
0137   //casj_final->SetBackgroundInputName("TowerBackground_Sub2");
0138   //casj_final->SetJetInputName("AntiKt_Tower_HIRecoSeedsRaw_r02");
0139   //casj_final->SetJetOutputName("AntiKt_Tower_r02_Sub1_V2");
0140   //casj_final->Verbosity(verbosity);
0141   //se->registerSubsystem(casj_final);
0142
0143   //cout << "Started Subtracting Towers" << endl;
0144
0145   SubtractTowers *st = new SubtractTowers();
0146   st->SetFlowModulation(G4HIJETS::do_flow);
0147   st->Verbosity(verbosity);
0148   se->registerSubsystem(st);
0149   /*
0150   NegativeTowerRemover *ntr = new NegativeTowerRemover();
0151   ntr->SetFlowModulation(G4HIJETS::do_flow);
0152   ntr->Verbosity(verbosity);
0153   se->registerSubsystem(ntr);
0154   */
0155
0156   //JetReco *towerjetreco = new JetReco();
0157   towerjetreco = new JetReco();
0158   towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWER_SUB1));
0159   towerjetreco->add_input(new TowerJetInput(Jet::HCALIN_TOWER_SUB1));
0160   towerjetreco->add_input(new TowerJetInput(Jet::HCALOUT_TOWER_SUB1));
0161   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.2, verbosity), "AntiKt_Tower_r02_Sub1");
0162   //towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.3, verbosity), "AntiKt_Tower_r03_Sub1");
0163   //towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.4, verbosity), "AntiKt_Tower_r04_Sub1");
0164   //towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.5, verbosity), "AntiKt_Tower_r05_Sub1");
0165
0166   towerjetreco->set_algo_node("ANTIKT");
0167   towerjetreco->set_input_node("TOWER");
0168   towerjetreco->Verbosity(verbosity);
0169   se->registerSubsystem(towerjetreco);
0170
0171   MyJetAnalysis *myJetAnalysis = new MyJetAnalysis("AntiKt_Tower_r02_Sub1", "AntiKt_Truth_r02", "myjet_ppPileup_R02_" + outname + ".root");
0172   myJetAnalysis->setPtRange(15, 200);
0173   myJetAnalysis->setEtaRange(-1.1, 1.1);
0174   myJetAnalysis->setMindR(0.2);
0175   //myJetAnalysis->doECCprint(false);
0176   se->registerSubsystem(myJetAnalysis);
0177
0178   //cout << "Finished Constructing the Jets" << endl;
0179
0180   //MyJetAnalysis *myJetAnalysis_R04 = new MyJetAnalysis("AntiKt_Tower_r04", "AntiKt_Truth_r04", "output/myjet_ppPileup_R04_" + outname + ".root");
0181   //myJetAnalysis_R04->setPtRange(15, 200);
0182   //myJetAnalysis_R04->setEtaRange(-1.1, 1.1);
0183   //myJetAnalysis_R04->setMindR(0.2);
0184   //myJetAnalysis_R04->doECCprint(false);
0185   //se->registerSubsystem(myJetAnalysis_R04);
0186
0187   string inputFile = "DST_TRUTH_JET_pythia8_Jet30_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-00000.root";
0188   string inputFile0 = "DST_CALO_CLUSTER_pythia8_Jet30_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-00000.root";
0189   string inputFile1 = "DST_TRUTH_pythia8_Jet30_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-00000.root";
0190   string inputFile2 = "DST_BBC_G4HIT_pythia8_Jet30_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-00000.root";
0191
0192   //cout << "Got Strings" << endl;
0193
0194   Fun4AllInputManager *intrue2 = new Fun4AllDstInputManager("DSTtruth2");
0195   intrue2->AddListFile(filelisttruthjet,1);
0196   se->registerInputManager(intrue2);
0197
0198   Fun4AllInputManager *intrue = new Fun4AllDstInputManager("DSTtruth");
0199   intrue->AddListFile(filelisttruth,1);
0200   se->registerInputManager(intrue);
0201
0202   Fun4AllInputManager *in = new Fun4AllDstInputManager("DSTcalo");
0203   in->AddListFile(filelistcalo,1);
0204   se->registerInputManager(in);
0205
0206   Fun4AllInputManager *in2 = new Fun4AllDstInputManager("DSThit");
0207   in2->AddListFile(filelisthits,1);
0208   se->registerInputManager(in2);
0209
0210   //se->skip(134);
0211
0212   //se->run(-1);
0213   //cout << "Ran" << endl;
0214   se->run(n_event);
0215   se->skip(n_skip);
0216
0217   se->End();
0218
0219   gSystem->Exit(0);
0220   return 0;
0221 }
0222 #endif