JetBkgdSub/macro/Fun4All_JetBkgd.C

0001 #pragma once
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003 #include <fun4all/SubsysReco.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/Fun4AllInputManager.h>
0006 #include <fun4all/Fun4AllDstInputManager.h>
0007
0008 #include <fun4all/Fun4AllDstOutputManager.h>
0009 #include <fun4all/Fun4AllOutputManager.h>
0010
0011 #include <phool/PHRandomSeed.h>
0012 #include <phool/recoConsts.h>
0013
0014 #include <g4centrality/PHG4CentralityReco.h>
0015
0016 #include <jetbackground/RetowerCEMC.h>
0017 #include <jetbackground/FastJetAlgoSub.h>
0018
0019 #include <jetbase/JetReco.h>
0020 #include <jetbase/TowerJetInput.h>
0021 #include <jetbase/FastJetAlgo.h>
0022
0023 #include <HIJetReco.C>
0024
0025 #include <jetbkgdsub/JetBkgdSub.h>
0026
0027 R__LOAD_LIBRARY(libfun4all.so)
0028 R__LOAD_LIBRARY(libg4jets.so)
0029 R__LOAD_LIBRARY(libjetbackground.so)
0030 R__LOAD_LIBRARY(libjetbase.so)
0031 R__LOAD_LIBRARY(libg4centrality.so)
0032 R__LOAD_LIBRARY(libg4dst.so)
0033 R__LOAD_LIBRARY(libJetBkgdSub.so)
0034
0035
0036
0037 #endif
0038
0039
0040 void Fun4All_JetBkgd(
0041   const char *filelisttruth = "dst_truth_jet.list",
0042   const char *filelistcalo = "dst_calo_cluster.list",
0043     const char *filelistglobal = "dst_global.list",
0044   const char * outputfile = "output.root",
0045   const double jet_parameter = 0.4
0046 )
0047 {
0048
0049   //-----------------------------------
0050   // Fun4All server initialization
0051   //-----------------------------------
0052
0053   // create fun4all server
0054   Fun4AllServer *se = Fun4AllServer::instance();
0055   int verbosity = 0;
0056   se->Verbosity(verbosity);
0057   recoConsts *rc = recoConsts::instance();
0058
0059   // centrality
0060   PHG4CentralityReco *cent = new PHG4CentralityReco();
0061   cent->Verbosity(0);
0062   cent->GetCalibrationParameters().ReadFromFile("centrality", "xml", 0, 0, string(getenv("CALIBRATIONROOT")) + string("/Centrality/"));
0063   se->registerSubsystem( cent );
0064
0065   // retower CEMC
0066   RetowerCEMC *rcemc = new RetowerCEMC();
0067   rcemc->Verbosity(verbosity);
0068   rcemc->set_towerinfo(true);
0069   se->registerSubsystem(rcemc);
0070
0071
0072   //-----------------------------------
0073   // Jet reco
0074   //-----------------------------------
0075   // Enable::HIJETS_TRUTH=false;
0076   // HIJetReco();
0077
0078   // tower jets
0079   // create jetreco and jettruth node names
0080   string rawname = "AntiKt_Tower_r0" + to_string(int(jet_parameter * 10));
0081   JetReco *towerjetreco = new JetReco();
0082   towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER));
0083   towerjetreco->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO));
0084   towerjetreco->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO));
0085   towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, jet_parameter), rawname);
0086   towerjetreco->set_algo_node("ANTIKT");
0087   towerjetreco->set_input_node("TOWER");
0088   towerjetreco->Verbosity(verbosity);
0089   se->registerSubsystem(towerjetreco);
0090
0091   // ==============
0092   // Jet Bkgd Sub
0093   // ==============
0094   double etamin = -1.1 + jet_parameter;
0095   double etamax = 1.1 - jet_parameter;
0096   JetBkgdSub *myJetTree = new JetBkgdSub(jet_parameter,outputfile);
0097   myJetTree->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER));
0098   myJetTree->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO));
0099   myJetTree->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO));
0100   myJetTree->doIterative(false);
0101   myJetTree->doAreaSub(true);
0102   myJetTree->doMultSub(true);
0103   myJetTree->doTruth(true);
0104   myJetTree->setMinRecoPt(5.0); // only sets range for reco jets
0105   myJetTree->setEtaRange(etamin, etamax);
0106   myJetTree->setPtRange(0, 100); // only sets range for truth jets
0107   myJetTree->Verbosity(verbosity);
0108   se->registerSubsystem(myJetTree);
0109
0110   //-----------------------------------
0111   // Input managers
0112   //-----------------------------------
0113
0114   Fun4AllInputManager *intrue = new Fun4AllDstInputManager("DSTtruth");
0115   intrue->AddListFile(filelisttruth,1);
0116   se->registerInputManager(intrue);
0117
0118   Fun4AllInputManager *in2 = new Fun4AllDstInputManager("DSTcalo");
0119   in2->AddListFile(filelistcalo,1);
0120   se->registerInputManager(in2);
0121
0122   Fun4AllInputManager *in3 = new Fun4AllDstInputManager("DSTglobal");
0123   in3->AddListFile(filelistglobal,1);
0124   se->registerInputManager(in3);
0125
0126   //-----------------------------------
0127   // Run the analysis
0128   //-----------------------------------
0129
0130   se->run(10);
0131   se->End();
0132
0133   gSystem->Exit(0);
0134   return 0;
0135
0136 }